Behavioural results for the TSP and 3SAT

project_folder = "/Volumes/GoogleDrive/My Drive/Melbourne/UNIMELB/Research/Complexity Project/ICx3_new/Code/behavAnalysis/ICx3_supp/"
code_folder = paste0(project_folder, "Code/")
data_folder = paste0(project_folder, "Data/")

setwd(code_folder)
knitr::opts_knit$set(root.dir = code_folder)

#Import functions
source(file.path(code_folder,"DescriptiveFunctions.R"))

# Participant Data
folderModels = "/Volumes/GoogleDrive/My Drive/Melbourne/UNIMELB/Research/Complexity Project/ICx3_new/Code/behavAnalysis/Output/"

folderOut_figures = paste0(project_folder,"Output/Figures")
folderOut_tables = paste0(project_folder,"Output/Tables")

# Features of the tasks needed for the analysis
SATTaskMaxTime=110
nVars_sat=5
SATTaskMaxClicks = 20

TSPTaskMaxTime = 40


# Brms parameters
chains_brms = 4
cores_brms = min(chains_brms,detectCores())
seed_brms = 111
iters_high = 4000

## Setting up the basics
library(ggplot2)
#library(stargazer)
library(knitr)
library(dplyr)
library(ggsignif)
library(pander)
library(plotly)
library(reshape2)
library(Hmisc)
library(readr)
library(brms)
library(parallel)
library(sjstats)
library(bayesplot)
library(bmlm)
library(DiagrammeR)
library(tidyr)

#From easystats:
library(parameters)
library(see)
library(bayestestR)
library(performance)

library(forcats)
library(ggpol)

# Functions

save_summarise_model = function(model, modelName){
  if (modelName %in% names(allModels)){
    warning("Model Name already exists!")
    model = allModels[[modelName]]
  } else {
    allModels[[modelName]]<<- model
  }
    table = model_parameters(model, 
                              ci_method = "HDI",ci=0.95,test = c("hdi","pd"))
    
    table2 = performance::model_performance(model,metrics="common")

    plo = plot(hdi(model, ci = c(0.95)), data = model)+
          ggtitle("Posterior distributions",
                  "with medians and 95% quantile-based-intervals") +
          theme(plot.title = element_text(hjust = 0.5),
                plot.subtitle = element_text(hjust = 0.5))
  
    print(as_tibble(table))
    print(as_tibble(table2))
    print(plo)
  
    #allModels[[modelName]]<<- model
}

# Use current models or rerun everything Comment one option

# allModels := list of all the regressions that are run 


## Option1: Use current saved models

# Import estimated models
allModels = read_rds(file.path(folderModels,"ICx3_models.RData"))

# Don't run brms models
brm2 =  function(...){
  return(FALSE)
}

# ## Option2: Rerun everything from scratch
# 
# allModels=vector("list", length=0)
# 
# brm2 =  function(...){
#   return(brm(...))
# }

1 Boolean satisfiability problem (3SAT)

1.1 Functions

calculate_sat_ICexpost =function(min_cost,nclauses){
  ICexpost  =  min_cost/nclauses
  return(ICexpost)
}

1.2 Accuracy

Read Clean Data to file

#Import SAT data
dataTrial_SAT= read_csv2(file.path(data_folder,"SAT_clean.csv"))

Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
  .default = col_double(),
  pID = col_character(),
  error = col_logical(),
  v = col_character(),
  l = col_character(),
  id = col_character(),
  region = col_character()
)
See spec(...) for full column specifications.

dataTrial_SAT_Time = read_csv2(file.path(data_folder,"SAT_clean_time.csv"))

Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
  .default = col_double(),
  pID = col_character(),
  error = col_logical(),
  v = col_character(),
  l = col_character(),
  id = col_character(),
  region = col_character()
)
See spec(...) for full column specifications.

dataTrial_SAT_Time$timeSpent_pct = dataTrial_SAT_Time$timeSpent/SATTaskMaxTime

color_scheme_set("green")

dataTrial_SAT$censored_time = ifelse(dataTrial_SAT$timeSpent == SATTaskMaxTime,
                             "right","none")

dataTrial_SAT_Time$censored_time = ifelse(dataTrial_SAT_Time$timeSpent == SATTaskMaxTime,
                             "right","none")

instanceProperties_SAT= read_rds(paste0(data_folder,"instance_properties_SAT.rds"))
instanceProperties_SAT = instanceProperties_SAT %>% select(instanceNumber, min_cost, min_model,costs,models,nSolutions)

dataTrial_SAT = left_join(dataTrial_SAT,instanceProperties_SAT, by = "instanceNumber")
dataTrial_SAT_Time = left_join(dataTrial_SAT_Time,instanceProperties_SAT, by = "instanceNumber")

dataTrial_SAT$ICexpost = calculate_sat_ICexpost(dataTrial_SAT$min_cost,dataTrial_SAT$nClauses)
dataTrial_SAT_Time$ICexpost = calculate_sat_ICexpost(dataTrial_SAT_Time$min_cost,dataTrial_SAT_Time$nClauses)

#CLicks
data_SAT_clicks = read_csv2(file.path(data_folder,"SAT_clicks.csv"))

Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
  pID = col_character(),
  block = col_double(),
  trial = col_double(),
  clicknumber = col_double(),
  Variable = col_double(),
  Literal = col_double(),
  state = col_double(),
  time = col_double()
)

nClicks = data_SAT_clicks %>% group_by(pID,block,trial) %>% summarise(n_clicks =n())

`summarise()` regrouping output by 'pID', 'block' (override with `.groups` argument)

dataTrial_SAT_clicks=merge(nClicks,
                  dataTrial_SAT, by=c("pID","block","trial"))

dataTrial_SAT_clicks$censored_clicks = ifelse(dataTrial_SAT_clicks$n_clicks== SATTaskMaxClicks,
                             "right","none")

#Summary Stats for Decision Problem
dataInput=dataTrial_SAT

accuracySummary = dataInput %>% group_by(pID) %>% summarise(acc=mean(correct)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() ungrouping output (override with .groups argument)

kable(accuracySummary, digits=2, caption = "Accuracy Summary")

Accuracy Summary

mean	min	max	SD
0.87	0.75	0.98	0.06

answerSummary = dataInput %>% group_by(pID) %>% summarise(acc=mean(answer)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() ungrouping output (override with .groups argument)

kable(answerSummary, digits=2, caption = "Proportion of times that YES was selected as the answer")

Proportion of times that YES was selected as the answer

mean	min	max	SD
0.45	0.28	0.61	0.09

yesNoProportions = dataInput %>% group_by(sol,pID) %>% summarise(acc=mean(correct)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() regrouping output by ‘sol’ (override with .groups argument) summarise() ungrouping output (override with .groups argument)

kable(yesNoProportions, digits=2, caption = "Accuracy By Solution")

Accuracy By Solution

sol	mean	min	max	SD
0	0.94	0.74	1	0.08
1	0.81	0.55	1	0.13

NA

1.2.1 Effect of Trial number on accuracy

#Trial (experience effect) effect
dataInput=dataTrial_SAT

#Summary
summaryByBlock = dataInput %>% group_by(block,pID) %>% summarise(acc=mean(correct)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

`summarise()` regrouping output by 'block' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

kable(summaryByBlock, digits=2, caption = "Accuracy By Block Number")

Accuracy By Block Number

block	mean	min	max	SD
1	0.85	0.75	1	0.07
2	0.87	0.69	1	0.08
3	0.88	0.69	1	0.09
4	0.87	0.69	1	0.10

#Regresssion
dataInput$totalTrial = dataInput$block * dataInput$trial

logitRandomIntercept = brm2(correct ~ totalTrial+ (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='sat_acc_01_r'

save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

1.2.2 TCC Accuracy Contrast

dataInput=dataTrial_SAT %>% group_by(phaseT,pID)%>%summarise(accuracy1=mean(correct))%>%ungroup()%>%group_by(phaseT)%>%summarise(accuracy=mean(accuracy1),se=se(accuracy1))%>%ungroup()

`summarise()` regrouping output by 'phaseT' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

dataInput$phaseT = recode(dataInput$phaseT, '0' = "Low IC", '1' = "High IC")

plo= ggplot(data=dataInput, aes(y=accuracy, x=as.factor(phaseT),label = round(accuracy,digits = 2))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1)+
  #geom_signif(comparisons = list(c("Low IC","High IC")), annotations="***", y_position = 0.95, tip_length = 0.03)+
  labs(title="Accuracy In and Out of phase Transition",x="Instance complexity",y="Accuracy")+
  coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white")+
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))

outputName = 'sat_acc_02_g'
plotExport(plo,outputName,folderOut_figures)

[1] "sat_acc_02_g"
quartz_off_screen 
                2 

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_SAT

#logistic with random effects (intercept)
logitRandomIntercept = brm2(correct ~ phaseT + (1|pID),
                           family=bernoulli(link="logit"),#binomial(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0)


tableName='sat_acc_02_r'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

1.2.2.1 Plot TCC

dataInput = dataTrial_SAT
dataInput = dataInput %>% 
  mutate(region = fct_relevel(region,
                              "Underconstrained", "Phase Transition", "Overconstrained"))%>% 
  mutate(sol= as.factor(sol)) %>% 
  mutate(sol = fct_relevel(sol,"1", "0"))

dataInput2 = dataInput %>% 
  group_by(instanceNumber,region,sol) %>% 
  summarise(accuracy=mean(correct))

`summarise()` regrouping output by 'instanceNumber', 'region' (override with `.groups` argument)

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7


plo = ggplot(dataInput2,aes(x=region,y=accuracy,fill=sol))+
  # geom_boxjitter(aes(fill=region),jitter.color = NA,jitter.shape = 21)+
  geom_boxjitter(jitter.color = NA,jitter.shape = 21,
                 jitter.params = list(height=0,seed=10),
                 outlier.shape= NA)+
                 #,outlier.shape = 4, outlier.size=0.9)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  scale_fill_manual(name  ="Region",
                    breaks=c("1","0"),
                    labels=c("Satisfiable","Unsatisfiable"),
                    values=c("#90BE6D","#F94144"))+
  xlab("Typical Case Complexity (TCC)")+
  ylab("Human Performance")+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(0.2,1)


plo
ggsave(paste0(folderOut_figures,"/TCC_SAT_acc.pdf"),plo,width = 7,height =6,units="in")

#brewer

1.2.3 Satisifability and TCC

dataInput=dataTrial_SAT %>% group_by(phaseT,sol,pID)%>%summarise(accuracy1=mean(correct))%>%ungroup()%>%group_by(sol,phaseT)%>%summarise(accuracy=mean(accuracy1),se=se(accuracy1))%>%ungroup()

`summarise()` regrouping output by 'phaseT', 'sol' (override with `.groups` argument)
`summarise()` regrouping output by 'sol' (override with `.groups` argument)

sol_names <- c(
                    `0` = "Correct answer: NO",
                    `1` = "Correct answer: YES",
                    `2` = "If this is here it means data not filtered"
                    )

plo=ggplot(data=dataInput, aes(y=accuracy, x=as.factor(as.logical(phaseT)), label = round(accuracy,digits = 2),group=1)) +
  geom_line()+
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1) +
  geom_point(shape=21, size=3, fill="white")+
  labs(title="Accuracy segregated by solvability",x="In Phase Transition?",y="Accuracy")+
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle=0,hjust=0.5))+
  geom_text(hjust=-0.5, colour="black")+
  facet_grid(.~sol, labeller = as_labeller(sol_names))

outputName = 'sat_acc_03_g'
plotExport(plo,outputName,folderOut_figures)

[1] "sat_acc_03_g"
quartz_off_screen 
                2 

dataInput= dataTrial_SAT
dataInput$phaseT = as.factor(dataInput$phaseT)
dataInput$sol = as.factor(dataInput$sol)
#logistic with random effects (intercept)
logitRandomIntercept = brm2(correct ~ phaseT + sol + phaseT:sol + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,                            seed=seed_brms,                            refresh = 0)

tableName='sat_acc_03_r'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

Marginal effects

#marginal_effects(logitRandomIntercept, ask=FALSE)
plot(marginal_effects(allModels[['sat_acc_03_r']]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

Effect of TCC on satisfiable instances

# st01
model = allModels[['sat_acc_03_r']]

fit = hypothesis(model,"phaseT1+ phaseT1:sol1=0", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |        95% HDI
--------------------------
H1        | [-2.53, -1.56]

1.2.3.1 Satisfiability alone

dataInput= dataTrial_SAT
dataInput$sol = as.factor(dataInput$sol)
#logistic with random effects (intercept)
logitRandomIntercept = brm2(correct ~ sol + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,                            seed=seed_brms,                            refresh = 0)

tableName='sat_acc_03B_r'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

1.2.4 Region effect on Accuracy

Overconstrained - Underconstrained - Satisfiability threshold

dataInput=dataTrial_SAT

# dataInput$pID=as.character(dataInput$pID)
# ps = sample(pull(dataInput,pID),20)
# dataInput = dataInput %>% filter(pID %in% ps)


dataInput2=dataInput %>% group_by(region,pID)%>%summarise(accuracy1=mean(correct))%>%
  ungroup()%>%group_by(region)%>%summarise(accuracy=mean(accuracy1),se=se(accuracy1))%>%ungroup()

`summarise()` regrouping output by 'region' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

dataInput2$region <- factor(dataInput2$region, levels = c('Underconstrained',
                                                          'Phase Transition',
                                                          'Overconstrained'))

plo= ggplot(data=dataInput2, aes(y=accuracy, x=region, label = round(accuracy,digits = 2))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1)+
  labs(title="Accuracy by region",x="Region",y="Accuracy")+
  coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white") +
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))

outputName = 'sat_acc_04_g'
plotExport(plo,outputName,folderOut_figures)

[1] "sat_acc_04_g"
quartz_off_screen 
                2 

#Stats test for time in and out of phase transition: ANOVA
dataInput$overConstrained= (dataInput$region=='Overconstrained')
dataInput$underConstrained= (dataInput$region=='Underconstrained')

logitRandomIntercept = brm2(correct ~ overConstrained + underConstrained + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='sat_acc_04_r_A'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

Is there a difference between accuracy in the overconstrained region and the underconstrained region?

fit = hypothesis(allModels[['sat_acc_04_r_A']],"overConstrainedTRUE=underConstrainedTRUE", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |       95% HDI
-------------------------
H1        | [-1.15, 0.24]

Region Effect on Accuracy by region and Satisfiability

dataInput=dataTrial_SAT 

# dataInput$pID=as.character(dataInput$pID)
# ps = sample(pull(dataInput,pID),20)
# dataInput = dataInput %>% filter(pID %in% ps)

dataInput2=dataInput %>% group_by(region,pID,sol)%>%summarise(accuracy1=mean(correct))%>%ungroup()%>%group_by(region,sol)%>%summarise(accuracy=mean(accuracy1),se=se(accuracy1))%>%ungroup()

`summarise()` regrouping output by 'region', 'pID' (override with `.groups` argument)
`summarise()` regrouping output by 'region' (override with `.groups` argument)

dataInput2$region <- factor(dataInput2$region, levels = c('Underconstrained',
                                                          'Phase Transition',
                                                          'Overconstrained'))

dataInput2$sol = recode(dataInput2$sol, '1' = "Solvable", '0' = "Non-Solvable")

plo= ggplot(data=dataInput2, aes(y=accuracy, x=region, label = round(accuracy,digits = 2))) +
  geom_bar(stat = "identity")+
  #geom_signif(annotations = c("***","***"),
  #            y_position=c(0.95,0.95),xmin=c(1.1,2.1),xmax=c(1.9,2.9))+
  #geom_signif(comparisons = list(c("Underconstrained","Overconstrained")),
  #             annotations="NS", y_position = 0.99, tip_length = 0.03)+
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1)+
  labs(title="Accuracy by region",x="Region",y="Accuracy")+
  coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white") +
  theme_light()+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))+
  facet_grid(.~sol)


outputName = 'sat_acc_04_g'
plotExport(plo,outputName,folderOut_figures)

[1] "sat_acc_04_g"
quartz_off_screen 
                2 

#Stats test for time in and out of phase transition: ANOVA
dataInput$overConstrained= (dataInput$region=='Overconstrained')
dataInput$underConstrained= (dataInput$region=='Underconstrained')

1.2.5 Number of solutions witnesses that satisfy the constraints

dataInput=dataTrial_SAT
#dataInput=nSolutions_sat(dataInput)

#Plots nSolutions (x) vs. Accuracy (y)
dataInput3 = dataInput %>% group_by(nSolutions,pID,phaseT)%>%summarise(accuracyMeans=mean(correct))%>%ungroup()%>%group_by(nSolutions,phaseT)%>%summarise(accuracy=mean(accuracyMeans),se=se(accuracyMeans))%>%ungroup()

`summarise()` regrouping output by 'nSolutions', 'pID' (override with `.groups` argument)
`summarise()` regrouping output by 'nSolutions' (override with `.groups` argument)

dataInput3$sol = dataInput3$nSolutions>=1
dataInput3$phaseT = recode(dataInput3$phaseT, '0' = "Low IC", '1' = "High IC")

plo= ggplot(data=dataInput3, aes(y=accuracy, x=as.factor(nSolutions))) +
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1)+
  geom_point(shape=23, size=3, fill="red")+
  labs(title="Accuracy and the Number of Solutions",
       x="Number of   Solutions",y="Accuracy")+
  coord_cartesian(ylim = c(0.5,1))+
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5),strip.text.x = element_blank())+
  facet_grid(phaseT~sol, scales = "free_x", space = "free_x")+
  geom_smooth(data=dataInput3, aes(y=accuracy, x=nSolutions),method=glm,se=FALSE,fullrange=TRUE,linetype = "dashed")

outputName = 'sat_acc_05_g'
plotExport(plo,outputName,folderOut_figures)

[1] "sat_acc_05_g"
quartz_off_screen 
                2 

Difference in the number of solutions of instances High/Low TCC (clearly, only for solvable instances)

dataInput2 = unique(dataInput %>% select(phaseT,id,nSolutions,sol) %>% filter(sol==1))
diffMeans = t.test(nSolutions ~ phaseT ,data=dataInput2)
pander(diffMeans)

-------------------------------------------------------------------
 Test statistic    df         P value       Alternative hypothesis 
---------------- ------- ----------------- ------------------------
     8.961        18.22   4.244e-08 * * *         two.sided        
-------------------------------------------------------------------

Table: Welch Two Sample t-test: `nSolutions` by `phaseT` (continued below)

 
-----------------------------------
 mean in group 0   mean in group 1 
----------------- -----------------
      6.588             1.312      
-----------------------------------

Regression for the number of witnesses analysis: correct ~ sol + phaseT + nSolutions + phaseT:nSolutions + (1|pID)

dataInput = dataInput %>% filter(sol==1)

#Includes a dummy if nSolutions==0, that's variable: sol.
logitRandomIntercept = brm2(correct ~ phaseT + nSolutions + phaseT:nSolutions + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,                            seed=seed_brms, refresh = 0)

tableName='sat_acc_05_r_A'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

This calculates the significance of the beta(slope) of the numberOfSolutions for high TCC instances.

#This version recodes phaseT to OutphaseT, to 

fit = hypothesis(allModels[['sat_acc_05_r_A']],"nSolutions + phaseT:nSolutions = 0", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |      95% HDI
------------------------
H1        | [0.64, 1.79]

Marginal Effects

plot(marginal_effects(allModels[['sat_acc_05_r_A']]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

1.2.5.1 NSol (regressed alone)

#st02
dataInput = dataTrial_SAT
dataInput = dataInput %>% filter(sol==1)

#Includes a dummy if nSolutions==0, that's variable: sol.
logitRandomIntercept = brm2(correct ~ nSolutions + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='sat_acc_05_r_B'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

mean_acc = dataInput %>% group_by(instanceNumber,nSolutions, sol,phaseT) %>%
  summarise(accuracy=mean(correct)) 

`summarise()` regrouping output by 'instanceNumber', 'nSolutions', 'sol' (override with `.groups` argument)

logitRandomIntercept = allModels[['sat_acc_05_r_B']]
pp  = plot(conditional_effects(logitRandomIntercept), plot = FALSE, ask = FALSE)

pp$nSolutions + 
  geom_point(data=mean_acc,aes(x = nSolutions, y = accuracy),inherit.aes = FALSE)

# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$nSolutions
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=mean_acc,
                               aes(x = nSolutions, y = accuracy, col=as.factor(sol),
                                   shape=as.factor(phaseT), size=2.5), inherit.aes = FALSE),
                         pp_plot$layers)


plo = pp_plot +
  xlab("Number of solution Witnesses")+#expression(IC[expost]))+
  ylab("Human Performance")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  scale_x_continuous(breaks = c(1,3,5,7,9),limits =c(1,10))+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(0.2,1)


plo
ggsave(paste0(folderOut_figures,"/Nwit_SAT_acc.pdf"),plo,width = 6,height =6,units="in")

1.2.5.2 NSol (with TCC and no interaction)

#st03
dataInput = dataTrial_SAT
dataInput = dataInput %>% filter(sol==1)

logitRandomIntercept = brm2(correct ~ nSolutions + phaseT + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='sat_acc_05_r_C'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

1.2.6 IC_expost

dataInput = dataTrial_SAT

# Summarise the data to plot
mean_accuracy = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(accuracy = mean(correct))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

ggplot(mean_accuracy, aes(x = ICexpost, y = accuracy))+
  geom_point() +
  theme_light() +
  stat_smooth(formula =  y ~ I(x^0.01), method="lm", se= FALSE)+
  xlab("IC_expost")+
  ylab("Human Accuracy")

dataInput = dataTrial_SAT
model_ICexpost = brm2(correct ~ ICexpost + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0)

tableName='sat_acc_icexpost_all'
save_summarise_model(model_ICexpost, tableName)

Model Name already exists!

mean_accuracy$correct= mean_accuracy$accuracy

model_ICexpost = allModels[['sat_acc_icexpost_all']]
pp  = plot(conditional_effects(model_ICexpost), plot = FALSE, ask = FALSE)

pp$ICexpost + 
  geom_point(data=mean_accuracy,aes(x = ICexpost, y = correct),inherit.aes = FALSE)

# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$ICexpost
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=mean_accuracy,
                               aes(x = ICexpost, y = correct, col=as.factor(sol),
                                   shape=as.factor(phaseT), size=2.5), inherit.aes = FALSE),
                         pp_plot$layers)


plo = pp_plot +
  xlab("Instance Complexity (IC)")+#expression(IC[expost]))+
  ylab("Human Performance")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(0.2,1)

plo
ggsave(paste0(folderOut_figures,"/IC_SAT_acc.pdf"),plo,width = 6,height =6,units="in")

# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

plo = pp$ICexpost +
geom_point(data=mean_accuracy,aes(x = ICexpost, y = correct, col=as.factor(phaseT),shape=as.factor(phaseT), size=2.5),inherit.aes = FALSE)+
  xlab("Instance Complexity (IC)")+#expression(IC[expost]))+
  ylab("Human Performance")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c( "#FC4E07","#E7B800"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))

plo
ggsave(paste0(folderOut_figures,"/IC_SAT_acc.pdf"),plo,width = 6,height =6,units="in")

Model Fit

model_ICexpost = allModels[['sat_acc_icexpost_all']]
dataInput = dataTrial_SAT

# Make predictions excluding random effects (pID)
predictions = predict(model_ICexpost,dataInput, re_formula = NA)


# Estimating the significance of the fit. This is done considering the probability estimation rather than the binary classification.

#Performs the Hosmer-Lemeshow goodness of fit test
logistic_significance = generalhoslem::logitgof(dataInput$correct, predictions[,1], g = 10, ord = FALSE)
logistic_significance

    Hosmer and Lemeshow test (binary model)

data:  dataInput$correct, predictions[, 1]
X-squared = 16.262, df = 8, p-value = 0.03878

# Finds R2 using binary outcomes
# https://stackoverflow.com/questions/40901445/function-to-calculate-r2-r-squared-in-r
#predictions = predict(model_ICexpost,dataInput, re_formula = NA)
#r_2_binary = cor(dataInput$ICexpost, predictions[,1])^2
r_2_binary = cor(dataInput$correct, predictions[,1])^2
r_2_binary

[1] 0.03659735

# Finds R2 using mean accuracies per instance
predictions2 = predict(model_ICexpost,mean_accuracy, re_formula = NA)
r_2_probabilities = cor(mean_accuracy$accuracy, predictions2[,1])^2
r_2_probabilities

[1] 0.1684351

1.2.6.1 Only unsatisfiable instances

dataInput = dataTrial_SAT
dataInput = dataInput %>% filter(sol==0)
model_ICexpost = brm2(correct ~ ICexpost + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0)

tableName='sat_acc_icexpost'
save_summarise_model(model_ICexpost, tableName)

Model Name already exists!

mean_accuracy2 = dataInput %>% 
  filter(sol == 0) %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(accuracy = mean(correct))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

mean_accuracy2$correct= mean_accuracy2$accuracy

model_ICexpost = allModels[['sat_acc_icexpost']]
pp  = plot(conditional_effects(model_ICexpost), plot = FALSE, ask = FALSE)

pp$ICexpost + 
  geom_point(data=mean_accuracy2,aes(x = ICexpost, y = correct),inherit.aes = FALSE)

Model Fit

model_ICexpost = allModels[['sat_acc_icexpost']]
#st05
dataInput = dataTrial_SAT %>% filter(sol==0)

# Make predictions excluding random effects (pID)
predictions = predict(model_ICexpost,dataInput, re_formula = NA)


# Estimating the significance of the fit. This is done considering the probability estimation rather than the binary classification.

#Performs the Hosmer-Lemeshow goodness of fit test
logistic_significance = generalhoslem::logitgof(dataInput$correct, predictions[,1], g = 10, ord = FALSE)
logistic_significance

    Hosmer and Lemeshow test (binary model)

data:  dataInput$correct, predictions[, 1]
X-squared = 32.849, df = 8, p-value = 6.557e-05

# Finds R2 using binary outcomes
# https://stackoverflow.com/questions/40901445/function-to-calculate-r2-r-squared-in-r
#predictions = predict(model_ICexpost,dataInput, re_formula = NA)
r_2_binary = cor(dataInput$correct, predictions[,1])^2
r_2_binary

[1] 0.0009535048

# Finds R2 using mean accuracies per instance
predictions2 = predict(model_ICexpost,mean_accuracy2, re_formula = NA)
r_2_probabilities = cor(mean_accuracy2$accuracy, predictions2[,1])^2
r_2_probabilities

[1] 0.001720168

1.3 Time Spent

The default data used in this section is dataTrial_SAT_Time. This dataset excludes those participants that never skipped to answer submission.

We first calculate some summary stats for the whole data set (dataTrial_SAT) and see how many participants were excluded in dataTrial_SAT_Time.

1.3.0.1 SummStats SAT

#Summary Stats for Decision Problem
dataInput=dataTrial_SAT

timeSummary = dataInput %>% group_by(pID) %>% summarise(acc=mean(timeSpent)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() ungrouping output (override with .groups argument)

kable(timeSummary, digits=1, caption = "Time Summary")

Time Summary

mean	min	max	SD
62.4	15.9	110	21.1

yesNoProportions = dataInput %>% group_by(sol,pID) %>% summarise(acc=mean(timeSpent)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() regrouping output by ‘sol’ (override with .groups argument) summarise() ungrouping output (override with .groups argument)

kable(yesNoProportions, digits=1, caption = "Time Spent By Solution")

Time Spent By Solution

sol	mean	min	max	SD
0	77.8	16.9	110	24.6
1	48.2	14.9	110	21.4

hist(dataInput$timeSpent,breaks=40)

all_pIDs = unique(dataTrial_SAT$pID)
pIDs_not_in_Time = all_pIDs[!(all_pIDs %in% unique(dataTrial_SAT_Time$pID))]
print(paste0("The participants removed in the TIME df are ", length(pIDs_not_in_Time)))

[1] "The participants removed in the TIME df are 1"

obs_removed = nrow(dataTrial_SAT)-nrow(dataTrial_SAT_Time)
print(paste0("The number of observations removed in the TIME df is ", obs_removed))

[1] "The number of observations removed in the TIME df is 63"

1.3.0.2 SummStats SAT_Time

#Summary Stats for Decision Problem
dataInput=dataTrial_SAT_Time

timeSummary = dataInput %>% group_by(pID) %>% summarise(acc=mean(timeSpent)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() ungrouping output (override with .groups argument)

kable(timeSummary, digits=1, caption = "Time Summary")

Time Summary

mean	min	max	SD
60.2	15.9	104.3	18.7

yesNoProportions = dataInput %>% group_by(sol,pID) %>% summarise(acc=mean(timeSpent)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() regrouping output by ‘sol’ (override with .groups argument) summarise() ungrouping output (override with .groups argument)

kable(yesNoProportions, digits=1, caption = "Time Spent By Solution")

Time Spent By Solution

sol	mean	min	max	SD
0	76.2	16.9	108.7	24.2
1	45.2	14.9	100.2	16.8

hist(dataInput$timeSpent,breaks=40)

1.3.0.3 Effect of Trial Number on Time-on-task

#Trial (experience effect) effect on timeSpent
dataInput=dataTrial_SAT_Time
dataInput$totalTrial = dataInput$block * dataInput$trial

summaryByBlock = dataInput %>% group_by(block,pID) %>% summarise(time=mean(timeSpent)) %>% summarise(mean=mean(time),min=min(time),max=max(time),SD=sd(time))

`summarise()` regrouping output by 'block' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

kable(summaryByBlock,digits=2 , caption="Time Spent per Trial segregated By Block")

Time Spent per Trial segregated By Block

block	mean	min	max	SD
1	71.95	19.38	110.00	21.27
2	63.90	16.81	110.00	19.56
3	55.52	14.69	110.00	22.26
4	49.34	1.68	87.13	21.73

#Regression
linearRandomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ totalTrial,
                            data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0, iter = iters_high,
                           control = list(adapt_delta = 0.98))

tableName='sat_time_01_r'
save_summarise_model(linearRandomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.

1.3.1 TCC

dataInput=dataTrial_SAT_Time %>% group_by(phaseT,pID)%>%summarise(timeSpent1=mean(timeSpent))%>%ungroup()%>%group_by(phaseT)%>%summarise(timeSpent=mean(timeSpent1),se=se(timeSpent1))%>%ungroup()

`summarise()` regrouping output by 'phaseT' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

dataInput$phaseT = recode(dataInput$phaseT, '0' = "Low IC", '1' = "High IC")

plo= ggplot(data=dataInput, aes(y=timeSpent, x=as.factor(phaseT),label = round(timeSpent,digits = 1))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=timeSpent-se, ymax=timeSpent+se), width=.1)+
  #geom_signif(comparisons = list(c("Low IC","High IC")), annotations="***", y_position = 0.95, tip_length = 0.03)+
  labs(title="Time Spent In and Out of phase Transition",x="Instance complexity",y="Time Spent")+
  #coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white")+
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))

outputName = "sat_time_02_g"
plotExport(plo,outputName,folderOut_figures)

[1] "sat_time_02_g"
quartz_off_screen 
                2 

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_SAT_Time

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(timeSpent| cens(censored_time) ~ phaseT + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_time_02_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

1.3.1.0.1 Plot TCC

dataInput = dataTrial_SAT_Time %>%
  mutate(region = fct_relevel(region,
                              "Underconstrained", "Phase Transition", "Overconstrained")) %>% 
  mutate(sol= as.factor(sol)) %>% 
  mutate(sol = fct_relevel(sol,"1", "0"))

dataInput$timeSpent_pct = dataInput$timeSpent/SATTaskMaxTime
dataInput2 = dataInput %>% 
  group_by(instanceNumber,region,sol, phaseT) %>% 
  summarise(timeSpent_pct=median(timeSpent_pct))

`summarise()` regrouping output by 'instanceNumber', 'region', 'sol' (override with `.groups` argument)

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

plo = ggplot(dataInput2,aes(x=region,y=timeSpent_pct,fill=sol))+
  geom_boxjitter(jitter.color = NA,jitter.shape = 21,
                 jitter.params = list(height=0,seed=10),
                 outlier.shape =NA)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  scale_fill_manual(name  ="Region",
                    breaks=c("1","0"),
                    labels=c("Satisfiable","Unsatisfiable"),
                    values=c("#90BE6D","#F94144"))+
  xlab("Typical Case Complexity (TCC)")+
  ylab("Time-on-task")+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  ylim(0.15,1)

plo
ggsave(paste0(folderOut_figures,"/TCC_SAT_time.pdf"),plo,width = 7,height =6,units="in")

1.3.2 Satsifiability and TCC

dataInput=dataTrial_SAT_Time %>% group_by(phaseT,sol,pID)%>%summarise(timeSpent1=mean(timeSpent))%>%ungroup()%>%group_by(sol,phaseT)%>%summarise(timeSpent=mean(timeSpent1),se=se(timeSpent1))%>%ungroup()

`summarise()` regrouping output by 'phaseT', 'sol' (override with `.groups` argument)
`summarise()` regrouping output by 'sol' (override with `.groups` argument)

sol_names <- c(
                    `0` = "Correct answer: NO",
                    `1` = "Correct answer: YES",
                    `2` = "If this is here it means data not filtered"
                    )

plo=ggplot(data=dataInput, aes(y=timeSpent, x=as.factor(as.logical(phaseT)), 
                               label = round(timeSpent,digits = 1),group=1)) +
  geom_line()+
  geom_errorbar(aes(ymin=timeSpent-se, ymax=timeSpent+se), width=.1) +
  geom_point(shape=21, size=3, fill="white")+
  labs(title="Time Spent segregated by solvability",x="In Phase Transition?",y="Time Spent")+
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle=0,hjust=0.5))+
  geom_text(hjust=-0.5, colour="black")+
  facet_grid(.~sol, labeller = as_labeller(sol_names))

outputName = "sat_time_03_g"
plotExport(plo,outputName,folderOut_figures)

[1] "sat_time_03_g"
quartz_off_screen 
                2 

dataInput= dataTrial_SAT_Time
#logistic with random effects (intercept)
dataInput$sol =as.factor(dataInput$sol)
dataInput$phaseT = as.factor(dataInput$phaseT)

# TTTODO
randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ phaseT + sol + phaseT:sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

title="Linear regressions"
tableName="sat_time_03_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

Marginal Effects

#marginal_effects(randomIntercept, ask=FALSE)
plot(marginal_effects(allModels[["sat_time_03_r"]]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

1.3.3 Satisfiability

#st06
dataInput= dataTrial_SAT_Time
#logistic with random effects (intercept)
#dataInput$sol =as.factor(dataInput$sol)

randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_time_03_r_B"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

1.3.4 Region effect on time-on-task (Constrainedness)

dataInput=dataTrial_SAT_Time #%>% mutate(region=recode(type, `1`= 'Phase Transition', `2`= 'Phase Transition', `3`= 'Phase Transition', `4`='Phase Transition', '5'='Overconstrained', '6'='Underconstrained'))

dataInput2=dataInput %>% group_by(region,pID)%>%summarise(timeSpent1=mean(timeSpent))%>%ungroup()%>%group_by(region)%>%summarise(timeSpent=mean(timeSpent1),se=se(timeSpent1))%>%ungroup()

`summarise()` regrouping output by 'region' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

dataInput2$region <- factor(dataInput2$region, levels = c('Underconstrained',
                                                          'Phase Transition',
                                                          'Overconstrained'))

plo= ggplot(data=dataInput2, aes(y=timeSpent, x=region, label = round(timeSpent,digits = 1))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=timeSpent-se, ymax=timeSpent+se), width=.1)+
  labs(title="time Spent by region",x="Region",y="time Spent")+
  #coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white") +
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))

outputName = "sat_time_04_g"
plotExport(plo,outputName,folderOut_figures)

[1] "sat_time_04_g"
quartz_off_screen 
                2 

#Stats test for time in and out of phase transition: ANOVA
dataInput$overConstrained= (dataInput$region=='Overconstrained')
dataInput$underConstrained= (dataInput$region=='Underconstrained')

randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ overConstrained + underConstrained + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_time_04_g_A"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

Is there a difference between time-on-task between the overconstrained region and the underconstrained region?

fit = hypothesis(allModels[["sat_time_04_g_A"]],"overConstrainedTRUE=underConstrainedTRUE", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |      95% HDI
------------------------
H1        | [0.38, 0.46]

1.3.5 Number of solution witnesses that satisfy the constraints

1.3.5.1 No. witnesses and TCC

dataInput=dataTrial_SAT_Time

#Plots nSolutions (x) vs. Accuracy (y)
dataInput3 = dataInput %>% group_by(nSolutions,pID,phaseT)%>%summarise(timeSpent1=mean(timeSpent))%>%ungroup()%>%group_by(nSolutions,phaseT)%>%summarise(timeSpent=mean(timeSpent1),se=se(timeSpent1))%>%ungroup()

`summarise()` regrouping output by 'nSolutions', 'pID' (override with `.groups` argument)
`summarise()` regrouping output by 'nSolutions' (override with `.groups` argument)

dataInput3$sol = dataInput3$nSolutions>=1
dataInput3$phaseT = recode(dataInput3$phaseT, '0' = "Low IC", '1' = "High IC")

plo= ggplot(data=dataInput3, aes(y=timeSpent, x=as.factor(nSolutions))) +
  geom_errorbar(aes(ymin=timeSpent-se, ymax=timeSpent+se), width=.1)+
  geom_point(shape=23, size=3, fill="red")+
  labs(title="Time Spent and the Number of Solutions",
       x="Number of   Solutions",y="Time Spent")+
  #coord_cartesian(ylim = c(0.5,1))+
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5),strip.text.x = element_blank())+
  facet_grid(phaseT~sol, scales = "free_x", space = "free_x")+
  geom_smooth(data=dataInput3, aes(y=timeSpent, x=nSolutions),method=glm,
              se=FALSE,fullrange=TRUE,linetype = "dashed")

outputName = "sat_time_05_g"
plotExport(plo,outputName,folderOut_figures)

[1] "sat_time_05_g"
quartz_off_screen 
                2 

# For satisfiable instance only
dataInput = dataInput %>% filter(sol==1)

randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ phaseT + nSolutions + phaseT:nSolutions + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_time_05_r_A"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

This calculates the significance of the beta(slope) of the number of witnesses for instances with high TCC.

#This version recodes phaseT to OutphaseT, to 

fit = hypothesis(allModels[["sat_time_05_r_A"]],"nSolutions + phaseT:nSolutions = 0", seed=seed_brms)

#fit = hypothesis(randomIntercept,"nSolutions = 0", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |        95% HDI
--------------------------
H1        | [-0.13, -0.04]

#summary(randomIntercept)

Marginal Effects

plot(marginal_effects(allModels[["sat_time_05_r_A"]]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

1.3.5.2 NSol (Regressed alone)

#st07
dataInput= dataTrial_SAT_Time
dataInput = dataInput %>% filter(sol==1)
#logistic with random effects (intercept)
#dataInput$sol =as.factor(dataInput$sol)

randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~  nSolutions + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_time_05_r_B"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

#dataInput3
mean_acc = dataInput %>% group_by(instanceNumber,nSolutions, sol,phaseT) %>%
  summarise(accuracy=mean(correct),
            timeSpent_med =median(timeSpent)) 

`summarise()` regrouping output by 'instanceNumber', 'nSolutions', 'sol' (override with `.groups` argument)

logitRandomIntercept = allModels[['sat_time_05_r_B']]
pp  = plot(conditional_effects(logitRandomIntercept,
           probs = c(0.025,0.975)), plot = FALSE, ask = FALSE)

pp$nSolutions + 
  geom_point(data=mean_acc,aes(x = nSolutions, y = timeSpent_med/SATTaskMaxTime),inherit.aes = FALSE)

NA
NA

# Improving the plot

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$nSolutions
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_boxplot(data=dataInput,
                                      aes(x = nSolutions, 
                                          group = nSolutions, 
                                          y = timeSpent_pct), 
                                      inherit.aes = FALSE),
                        geom_point(data=mean_acc, 
                                   aes(x = nSolutions, y = timeSpent_med/SATTaskMaxTime,
                                       col=as.factor(sol), shape=as.factor(phaseT), size=2.5),
                                   inherit.aes = FALSE),
                         pp_plot$layers)

plo = pp_plot +
  xlab("Number of solution Witnesses")+#expression(IC[expost]))+
  ylab("Time-on-task")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  ylim(0, 1)+
  scale_x_continuous(breaks = c(1,3,5,7,9),limits =c(0.5,10.5))+
  guides(shape = guide_legend(override.aes = list(size = 3)))

plo
ggsave(paste0(folderOut_figures,"/Nwit_SAT_time.pdf"),plo,width = 6,height =6,units="in")

# Improving the plot

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$nSolutions
pp_plot$layers[[1]]$geom_params$se = TRUE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=dataInput,
                               aes(x = nSolutions, y = timeSpent_pct,
                                   shape=as.factor(phaseT), size=0.4,
                                   alpha = 0.5),
                                   inherit.aes = FALSE),
                    geom_point(data=mean_acc, 
                                   aes(x = nSolutions, 
                                       y = timeSpent_med/SATTaskMaxTime,
                                       col=as.factor(sol),
                                       shape=as.factor(phaseT), size=2.5),
                                   inherit.aes = FALSE),
                         pp_plot$layers)


plo = pp_plot +
  xlab("Number of solution Witnesses")+#expression(IC[expost]))+
  ylab("Time-on-task")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+#c( "#FC4E07","#E7B800")
  # scale_shape_manual(name="IC",values = c(2, 8))+
  # scale_color_manual(name="Solution",values = c("red", "blue"))+
  theme_classic()+

  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  ylim(0, 1)+
  scale_x_continuous(breaks = c(1,3,5,7,9),limits =c(0.5,10.5))+
  guides(shape = guide_legend(override.aes = list(size = 3)))

plo
ggsave(paste0(folderOut_figures,"/Nwit_SAT_time2.pdf"),plo,width = 6,height =6,units="in")

1.3.6 IC_expost (distance from the optimum)

dataInput = dataTrial_SAT_Time
# Summarise the data to plot
mean_accuracy_time = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(accuracy = mean(correct), 
            timeSpent_med = median(timeSpent),
            timeSpent = mean(timeSpent))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

ggplot(mean_accuracy_time,aes(x = ICexpost, y = timeSpent))+
  geom_point() +
  theme_light() +
  stat_smooth(formula =  y ~ I(x^0.01), method="lm", se= FALSE)+
  xlab("IC_expost")+
  ylab("Time Spent")

1.3.6.0.1 Unsatisfiable instances

dataInput = dataTrial_SAT_Time %>% filter(sol==0)

model_ICexpost = brm2(timeSpent_pct | cens(censored_time) ~ ICexpost + (1|pID),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0, iter = iters_high)

tableName='sat_time_icexpost'
save_summarise_model(model_ICexpost, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

1.3.6.0.2 All instances

dataInput = dataTrial_SAT_Time

model_ICexpost = brm2(timeSpent_pct | cens(censored_time) ~ ICexpost + (1|pID),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0, iter = iters_high)

tableName='sat_time_icexpost2'
save_summarise_model(model_ICexpost, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

Marginal Effects

mean_accuracy_time$correct= mean_accuracy_time$accuracy
#mean_accuracy_time$ICexpost05 = mean_accuracy_time$ICexpost^0.5

model_ICexpost = allModels[['sat_time_icexpost2']]
pp  = plot(conditional_effects(model_ICexpost), plot = FALSE, ask = FALSE)

pp$ICexpost + 
  geom_point(data=mean_accuracy_time,aes(x = ICexpost, y = timeSpent/SATTaskMaxTime),inherit.aes = FALSE)

# Improving the plot option 2

pp_plot = pp$ICexpost
pp_plot$layers[[1]]$geom_params$se = TRUE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=dataInput,
                               aes(x = ICexpost, y = timeSpent_pct,
                                   shape=as.factor(phaseT), size=0.4,
#                                   col=as.factor(sol),
                                   alpha = 0.2),
                                   inherit.aes = FALSE),
                    geom_point(data=mean_accuracy_time, 
                                   aes(x = ICexpost, 
                                       y = timeSpent_med/SATTaskMaxTime,
                                       col=as.factor(sol),
                                       shape=as.factor(phaseT), size=2.5),
                                   inherit.aes = FALSE),
                         pp_plot$layers)

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

plo = pp_plot +
  xlab("Instance Complexity (IC)")+#expression(IC[expost]))+
  ylab("Time-on-task")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+#c( "#FC4E07","#E7B800")
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  ylim(0, 1.1)+
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1))+
  guides(shape = guide_legend(override.aes = list(size = 3)))

Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.

plo
ggsave(paste0(folderOut_figures,"/IC_SAT_time1.pdf"),plo,width = 6,height =6,units="in")

1.4 Time-on-task and Accuracy

dataInput= dataTrial_SAT_Time %>% group_by(correct,pID) %>% summarise(timeSpentMeans=mean(timeSpent)) %>% ungroup() %>% group_by(correct) %>%
  summarise(time=mean(timeSpentMeans),se=se(timeSpentMeans))%>%ungroup()

`summarise()` regrouping output by 'correct' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

plo = ggplot(data=dataInput, aes(y=time, x=as.factor(as.logical(correct)), label = round(time,digits = 2),group=1)) +
  geom_line()+
  geom_errorbar(aes(ymin=time-se, ymax=time+se), width=.1) +
  geom_point(shape=21, size=3, fill="white")+
  labs(title="Time Spent on Correct/Incorrect instances",x="Answer Correct?",y="Time Spent") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle=0,hjust=0.5))+
  geom_text(hjust=-0.5, colour="black")

outputName = "sat_time_09_g"
plotExport(plo,outputName,folderOut_figures)

[1] "sat_time_09_g"
quartz_off_screen 
                2 

dataInput = dataTrial_SAT_Time
#Stats test for time for correct/incorrect answers
logitRandomIntercept = brm2(correct ~ timeSpent_pct + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)


tableName="sat_time_09_r"
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

1.5 Number of Clicks

1.5.0.1 Satisfiability

dataInput=dataTrial_SAT_clicks

# BOXPLOT
plo= ggplot(data=dataInput, aes(y=n_clicks, x=region, label = round(n_clicks,digits = 0))) +
  geom_boxplot()+
  labs(x="Region",y="n_clicks")+
  #coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white") +
  theme_light()+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))+
  facet_grid(.~sol)


outputName = 'sat_acc_04_g'
plotExport(plo,outputName,folderOut_figures)

[1] "sat_acc_04_g"
quartz_off_screen 
                2 

Nice Plot:

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

dataInput2 = dataTrial_SAT_clicks %>% 
  mutate(region = fct_relevel(region,
                              "Underconstrained", "Phase Transition", "Overconstrained"))%>% 
  mutate(sol= as.factor(sol)) %>% 
  mutate(sol = fct_relevel(sol,"1", "0"))

dataInput3 = dataInput2 %>% 
  group_by(instanceNumber,region,sol) %>% 
  summarise(n_clicks=mean(n_clicks))

`summarise()` regrouping output by 'instanceNumber', 'region' (override with `.groups` argument)

plo = ggplot(dataInput3,aes(x=region,y=n_clicks,fill=sol))+
  # geom_boxjitter(aes(fill=region),jitter.color = NA,jitter.shape = 21)+
  geom_boxjitter(jitter.color = NA,jitter.shape = 21,
                 jitter.params = list(height=0,seed=10),
                 outlier.shape= NA)+
                 #,outlier.shape = 4, outlier.size=0.9)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  scale_fill_manual(name  ="Region",
                    breaks=c("1","0"),
                    labels=c("Satisfiable","Unsatisfiable"),
                    values=c("#90BE6D","#F94144"))+
  xlab("Typical Case Complexity (TCC)")+
  ylab("Number of clicks")+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(0,20)


plo
ggsave(paste0(folderOut_figures,"/TCC_SAT_nclicks.pdf"),plo,width = 7,height =6,units="in")

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_SAT_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks| cens(censored_clicks) ~ sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_clicks_01_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

1.5.0.2 TCC

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_SAT_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks| cens(censored_clicks) ~ phaseT + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_clicks_02_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

1.5.0.2.1 And Satisfiability

dataInput= dataTrial_SAT_clicks

#logistic with random effects (intercept)
dataInput$phaseT = as.factor(dataInput$phaseT)
dataInput$sol = as.factor(dataInput$sol)
randomIntercept = brm2(n_clicks| cens(censored_clicks) ~ phaseT + sol + phaseT:sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_clicks_02B_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

pp  = plot(conditional_effects(allModels[["sat_clicks_02B_r"]]), plot = TRUE, ask = FALSE)

1.5.0.3 ICexpost

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_SAT_clicks

#logistic with random effects (intercept)
randomIntercept = brm2(n_clicks| cens(censored_clicks) ~ ICexpost + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_clicks_03_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

dataInput = dataTrial_SAT_clicks
# Summarise the data to plot
mean_nclicks = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(n_clicks = median(n_clicks))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

model_ICexpost = allModels[["sat_clicks_03_r"]]
pp  = plot(conditional_effects(model_ICexpost), plot = FALSE, ask = FALSE)

pp$ICexpost +
  geom_point(data=mean_nclicks,aes(x = ICexpost, y =n_clicks),inherit.aes = FALSE)

# Improving the plot
# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$ICexpost
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=mean_nclicks,
                               aes(x = ICexpost, y = n_clicks, col=as.factor(sol),
                                   shape=as.factor(phaseT), size=2.5), inherit.aes = FALSE),
                         pp_plot$layers)


plo = pp_plot +
  xlab("Instance Complexity (IC)")+#expression(IC[expost]))+
  ylab("Number of Clicks")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))

plo
ggsave(paste0(folderOut_figures,"/IC_SAT_nclicks.pdf"),plo,width = 6,height =6,units="in")

Only Unsatisfiable:

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_SAT_clicks %>% filter(sol==0)

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks| cens(censored_clicks) ~ ICexpost + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_clicks_03_rB"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

dataInput= dataTrial_SAT_clicks #%>% filter(sol==0)
# Summarise the data to plot
mean_nclicks = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(n_clicks = median(n_clicks))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

model_ICexpost = allModels[["sat_clicks_03_rB"]]
pp  = plot(conditional_effects(model_ICexpost), plot = FALSE, ask = FALSE)

pp$ICexpost +
  geom_point(data=mean_nclicks,aes(x = ICexpost, y =n_clicks),inherit.aes = FALSE)

# Improving the plot
# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$ICexpost
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_jitter(data=dataInput,
                               aes(x = ICexpost, y = n_clicks,
                                   shape=as.factor(phaseT), size=0.4,
#                                   col=as.factor(sol),
                                   alpha = 0.2),
                                   inherit.aes = FALSE,height=0.2),
                        geom_point(data=mean_nclicks, 
                                   aes(x = ICexpost, y = n_clicks,
                                       col=as.factor(sol), shape=as.factor(phaseT), size=2.5),
                                   inherit.aes = FALSE),
                         pp_plot$layers)

plo = pp_plot +
  xlab("Instance Complexity (IC)")+#expression(IC[expost]))+
  ylab("Number of Clicks")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+ 
  xlim(0,0.09)

plo
ggsave(paste0(folderOut_figures,"/IC_SAT_nclicks_unsat.pdf"),plo,width = 6,height =6,units="in")

1.5.0.3.1 And Satisfiability

dataInput= dataTrial_SAT_clicks

#logistic with random effects (intercept)
dataInput$unsol = abs(dataInput$sol-1)
randomIntercept = brm2(n_clicks| cens(censored_clicks) ~ unsol + unsol:ICexpost + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_clicks_03_rC"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

pp  = plot(conditional_effects(allModels[["sat_clicks_03_rC"]]), plot = TRUE, ask = FALSE)

1.5.0.4 Nsolutions

dataInput= dataTrial_SAT_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks| cens(censored_clicks) ~ nSolutions + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_clicks_04_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

Only Satisfiable:

dataInput= dataTrial_SAT_clicks %>% filter(sol==1)

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks| cens(censored_clicks) ~ nSolutions + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="sat_clicks_04_rB"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

#dataInput = dataTrial_SAT_clicks
# Summarise the data to plot
mean_nclicks = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(n_clicks = median(n_clicks))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

model_ICexpost = allModels[["sat_clicks_04_rB"]]
pp  = plot(conditional_effects(model_ICexpost), plot = FALSE, ask = FALSE)

pp$nSolutions +
  geom_point(data=mean_nclicks,aes(x = nSolutions, y =n_clicks),inherit.aes = FALSE)

# Improving the plot
# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$nSolutions
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_jitter(data=dataInput,
                               aes(x = nSolutions, y = n_clicks,
                                   shape=as.factor(phaseT), size=0.4,
#                                   col=as.factor(sol),
                                   alpha = 0.2),
                                   inherit.aes = FALSE,width=0.2,height =0.2),
                        geom_point(data=mean_nclicks, 
                                   aes(x = nSolutions, y = n_clicks,
                                       col=as.factor(sol), shape=as.factor(phaseT), size=2.5),
                                   inherit.aes = FALSE),
                         pp_plot$layers)

plo = pp_plot +
  xlab("Number of Solution Witnesses")+
  ylab("Number of Clicks")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))

plo
ggsave(paste0(folderOut_figures,"/Nwit_SAT_nclicks.pdf"),plo,width = 6,height =6,units="in")

2 Traveling salesperson Problem

2.1 Functions

calculate_tsp_ICexpost =function(opt,max_dist){
  ICexpost  =  abs(opt- max_dist)/opt
  return(ICexpost)
}

2.2 Accuracy

Read Clean Data to file

dataTrial_TSP= read_csv2(file.path(data_folder,"TSP_clean.csv"))

Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
  .default = col_double(),
  pID = col_character(),
  itemsSelected = col_character(),
  error = col_logical(),
  cx = col_character(),
  cy = col_character(),
  distances = col_character(),
  id = col_character(),
  region = col_character()
)
See spec(...) for full column specifications.

dataTrial_TSP_Time = read_csv2(file.path(data_folder,"TSP_clean_time.csv"))

Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
  .default = col_double(),
  pID = col_character(),
  itemsSelected = col_character(),
  error = col_logical(),
  cx = col_character(),
  cy = col_character(),
  distances = col_character(),
  id = col_character(),
  region = col_character()
)
See spec(...) for full column specifications.

dataTrial_TSP_Time$timeSpent_pct = dataTrial_TSP_Time$timeSpent/TSPTaskMaxTime

color_scheme_set("purple")

dataTrial_TSP$censored_time = ifelse(dataTrial_TSP$timeSpent == TSPTaskMaxTime,
                             "right","none")

dataTrial_TSP_Time$censored_time = ifelse(dataTrial_TSP_Time$timeSpent == TSPTaskMaxTime,
                             "right","none")

instanceProperties_TSP= read_csv2(paste0(data_folder,"instance_properties_TSP.csv"))

Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
  instanceNumber = col_double(),
  nSolutions = col_double(),
  opt_distance = col_double(),
  opt_tour = col_character(),
  cx = col_character(),
  cy = col_character(),
  distances = col_character(),
  id = col_character(),
  type = col_double(),
  sol = col_double(),
  max_distance = col_double(),
  nCities = col_double(),
  param = col_character(),
  dist_from_opt = col_double()
)

instanceProperties_TSP =instanceProperties_TSP %>% 
  select(instanceNumber,nSolutions,opt_distance,opt_tour,dist_from_opt)

dataTrial_TSP = left_join(dataTrial_TSP,instanceProperties_TSP, by = "instanceNumber")
dataTrial_TSP_Time = left_join(dataTrial_TSP_Time,instanceProperties_TSP, by = "instanceNumber")

dataTrial_TSP$ICexpost = calculate_tsp_ICexpost(dataTrial_TSP$opt_distance,
                                          dataTrial_TSP$max_distance)

dataTrial_TSP_Time$ICexpost = calculate_tsp_ICexpost(dataTrial_TSP_Time$opt_distance,
                                          dataTrial_TSP_Time$max_distance)


data_TSP_clicks = read_csv2(file.path(data_folder,"TSP_clicks.csv"))

Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
  pID = col_character(),
  block = col_double(),
  trial = col_double(),
  citynumber.100.Reset. = col_double(),
  In.1..Out.0..Reset.3. = col_double(),
  time = col_double()
)

nClicks = data_TSP_clicks %>% group_by(pID,block,trial) %>% summarise(n_clicks =n())

`summarise()` regrouping output by 'pID', 'block' (override with `.groups` argument)

#View(nClicks)

dataTrial_TSP_clicks=merge(nClicks,
                  dataTrial_TSP, by=c("pID","block","trial"))

#Summary Stats for Decision Problem
dataInput=dataTrial_TSP

accuracySummary = dataInput %>% group_by(pID) %>% summarise(acc=mean(correct)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() ungrouping output (override with .groups argument)

kable(accuracySummary, digits=2, caption = "Accuracy Summary")

Accuracy Summary

mean	min	max	SD
0.85	0.76	0.93	0.05

answerSummary = dataInput %>% group_by(pID) %>% summarise(acc=mean(answer)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() ungrouping output (override with .groups argument)

kable(answerSummary, digits=2, caption = "Proportion of times that YES was selected as the answer")

Proportion of times that YES was selected as the answer

mean	min	max	SD
0.5	0.35	0.68	0.09

yesNoProportions = dataInput %>% group_by(sol,pID) %>% summarise(acc=mean(correct)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() regrouping output by ‘sol’ (override with .groups argument) summarise() ungrouping output (override with .groups argument)

kable(yesNoProportions, digits=2, caption = "Accuracy By Solution")

Accuracy By Solution

sol	mean	min	max	SD
0	0.86	0.63	1	0.11
1	0.85	0.69	1	0.09

NA

2.2.1 Effect of trial number on accuracy

#Trial (experience effect) effect
dataInput=dataTrial_TSP

#Summary
summaryByBlock = dataInput %>% group_by(block,pID) %>% summarise(acc=mean(correct)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

`summarise()` regrouping output by 'block' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

kable(summaryByBlock, digits=2, caption = "Accuracy By Block Number")

Accuracy By Block Number

block	mean	min	max	SD
1	0.83	0.58	0.96	0.09
2	0.86	0.75	0.96	0.06
3	0.87	0.71	1.00	0.08

#Regresssion
dataInput$totalTrial = dataInput$block * dataInput$trial

logitRandomIntercept = brm2(correct ~ totalTrial+ (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0)

tableName='tsp_acc_01_r'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

2.2.2 TCC

dataInput=dataTrial_TSP %>% group_by(phaseT,pID)%>%summarise(accuracy1=mean(correct))%>%ungroup()%>%group_by(phaseT)%>%summarise(accuracy=mean(accuracy1),se=se(accuracy1))%>%ungroup()

`summarise()` regrouping output by 'phaseT' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

dataInput$phaseT = recode(dataInput$phaseT, '0' = "Low IC", '1' = "High IC")

plo= ggplot(data=dataInput, aes(y=accuracy, x=as.factor(phaseT),label = round(accuracy,digits = 2))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1)+
  #geom_signif(comparisons = list(c("Low IC","High IC")), annotations="***", y_position = 0.95, tip_length = 0.03)+
  labs(title="Accuracy In and Out of phase Transition",x="Instance complexity",y="Accuracy")+
  coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white")+
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))

outputName = 'tsp_acc_02_g'
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_acc_02_g"
quartz_off_screen 
                2 

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP

#logistic with random effects (intercept)
logitRandomIntercept = brm2(correct ~ phaseT + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='tsp_acc_02_r'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

2.2.2.1 Plot TCC

dataInput = dataTrial_TSP

dataInput = dataInput %>% 
  mutate(region = fct_relevel(region,
                              "Underconstrained", "Phase Transition", "Overconstrained")) %>% 
  mutate(sol= as.factor(sol)) %>% 
  mutate(sol = fct_relevel(sol,"1", "0"))

dataInput2 = dataInput %>% 
  group_by(instanceNumber,region,sol) %>% 
  summarise(accuracy=mean(correct))

`summarise()` regrouping output by 'instanceNumber', 'region' (override with `.groups` argument)

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7


plo = ggplot(dataInput2,aes(x=region,y=accuracy,fill=sol))+
  # geom_boxjitter(aes(fill=region),jitter.color = NA,jitter.shape = 21)+
  geom_boxjitter(jitter.color = NA,jitter.shape = 21,
                 jitter.params = list(height=0,seed=10),
                 outlier.shape= NA)+
                 #,outlier.shape = 4, outlier.size=0.9)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  scale_fill_manual(name  ="Region",
                    breaks=c("1","0"),
                    labels=c("Satisfiable","Unsatisfiable"),
                    values=c("#90BE6D","#F94144"))+
  xlab("Typical Case Complexity (TCC)")+
  ylab("Human Performance")+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(0.2,1)


plo
ggsave(paste0(folderOut_figures,"/TCC_TSP_acc.pdf"),plo,width = 7,height =6,units="in")

2.2.3 Satsifiability and TCC

dataInput=dataTrial_TSP %>% group_by(phaseT,sol,pID)%>%summarise(accuracy1=mean(correct))%>%ungroup()%>%group_by(sol,phaseT)%>%summarise(accuracy=mean(accuracy1),se=se(accuracy1))%>%ungroup()

`summarise()` regrouping output by 'phaseT', 'sol' (override with `.groups` argument)
`summarise()` regrouping output by 'sol' (override with `.groups` argument)

sol_names <- c(
                    `0` = "Correct answer: NO",
                    `1` = "Correct answer: YES",
                    `2` = "If this is here it means data not filtered"
                    )

plo=ggplot(data=dataInput, aes(y=accuracy, x=as.factor(as.logical(phaseT)), label = round(accuracy,digits = 2),group=1)) +
  geom_line()+
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1) +
  geom_point(shape=21, size=3, fill="white")+
  labs(title="Accuracy segregated by solvability",x="In Phase Transition?",y="Accuracy")+
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle=0,hjust=0.5))+
  geom_text(hjust=-0.5, colour="black")+
  facet_grid(.~sol, labeller = as_labeller(sol_names))

outputName = 'tsp_acc_03_g'
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_acc_03_g"
quartz_off_screen 
                2 

dataInput= dataTrial_TSP

dataInput$phaseT = as.factor(dataInput$phaseT)
dataInput$sol = as.factor(dataInput$sol)

#logistic with random effects (intercept)
logitRandomIntercept = brm2(correct ~ phaseT + sol + phaseT:sol + (1|pID), 
                           family=bernoulli(link="logit"),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='tsp_acc_03_r'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

Marginal effects

plot(marginal_effects(allModels[['tsp_acc_03_r']]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

Effect of TCC on satisfiable instances

#st01
model = allModels[['tsp_acc_03_r']]

fit = hypothesis(model,"phaseT1+ phaseT1:sol1=0", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |        95% HDI
--------------------------
H1        | [-2.61, -1.52]

2.2.3.1 Only satisfiability

dataInput= dataTrial_TSP
dataInput$sol = as.factor(dataInput$sol)

#logistic with random effects (intercept)
logitRandomIntercept = brm2(correct ~ sol + (1|pID), 
                           family=bernoulli(link="logit"),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='tsp_acc_03B_r'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

2.2.4 Region effect on Accuracy (Constrainedness)

dataInput=dataTrial_TSP

dataInput2=dataInput %>% group_by(region,pID)%>%summarise(accuracy1=mean(correct))%>%ungroup()%>%group_by(region)%>%summarise(accuracy=mean(accuracy1),se=se(accuracy1))%>%ungroup()

`summarise()` regrouping output by 'region' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

dataInput2$region <- factor(dataInput2$region, levels = c('Underconstrained',
                                                          'Phase Transition',
                                                          'Overconstrained'))

plo= ggplot(data=dataInput2, aes(y=accuracy, x=region, label = round(accuracy,digits = 2))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1)+
  labs(title="Accuracy by region",x="Region",y="Accuracy")+
  coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white") +
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))

outputName = 'tsp_acc_04_g'
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_acc_04_g"
quartz_off_screen 
                2 

#Stats test for time in and out of phase transition: ANOVA
dataInput$overConstrained= (dataInput$region=='Overconstrained')
dataInput$underConstrained= (dataInput$region=='Underconstrained')

logitRandomIntercept = brm2(correct ~ overConstrained + underConstrained + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,                            seed=seed_brms, refresh = 0)

tableName='tsp_acc_04_r_A'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

Is there a difference between accuracy in the overconstrained region and the underconstrained region?

fit = hypothesis(allModels[['tsp_acc_04_r_A']],"overConstrainedTRUE=underConstrainedTRUE", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |       95% HDI
-------------------------
H1        | [-0.61, 0.84]

Region effect on accuracy by region and satisfiability

dataInput=dataTrial_TSP 

# dataInput$pID=as.character(dataInput$pID)
# ps = sample(pull(dataInput,pID),20)
# dataInput = dataInput %>% filter(pID %in% ps)

dataInput2=dataInput %>% group_by(region,pID,sol)%>%summarise(accuracy1=mean(correct))%>%ungroup()%>%group_by(region,sol)%>%summarise(accuracy=mean(accuracy1),se=se(accuracy1))%>%ungroup()

`summarise()` regrouping output by 'region', 'pID' (override with `.groups` argument)
`summarise()` regrouping output by 'region' (override with `.groups` argument)

dataInput2$region <- factor(dataInput2$region, levels = c('Underconstrained',
                                                          'Phase Transition',
                                                          'Overconstrained'))

dataInput2$sol = recode(dataInput2$sol, '1' = "Solvable", '0' = "Non-Solvable")

plo= ggplot(data=dataInput2, aes(y=accuracy, x=region, label = round(accuracy,digits = 2))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1)+
  labs(title="Accuracy by region",x="Region",y="Accuracy")+
  coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white") +
  theme_light()+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))+
  facet_grid(.~sol)


outputName = 'tsp_acc_05_g'
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_acc_05_g"
quartz_off_screen 
                2 

#Stats test for time in and out of phase transition: ANOVA
dataInput$overConstrained= (dataInput$region=='Overconstrained')
dataInput$underConstrained= (dataInput$region=='Underconstrained')

2.2.5 Number of solution witnesses

Difference in the number of witnesses for instances with low/high TCC

#Number of solutions of solvable solutions: Mean difference between low/high TCC

dataInput2 = unique(dataInput %>% select(phaseT,id,nSolutions,sol) %>% filter(sol==1))
diffMeans = t.test(nSolutions ~ phaseT ,data=dataInput2)
pander(diffMeans)

-------------------------------------------------------------------
 Test statistic    df         P value       Alternative hypothesis 
---------------- ------- ----------------- ------------------------
     7.264        32.84   2.543e-08 * * *         two.sided        
-------------------------------------------------------------------

Table: Welch Two Sample t-test: `nSolutions` by `phaseT` (continued below)

 
-----------------------------------
 mean in group 0   mean in group 1 
----------------- -----------------
      26208             4720       
-----------------------------------

dataInput=dataTrial_TSP

#Plots nSolutions (x) vs. Accuracy (y)
dataInput3 = dataInput %>% group_by(nSolutions,pID,phaseT)%>%summarise(accuracyMeans=mean(correct))%>%ungroup()%>%group_by(nSolutions,phaseT)%>%summarise(accuracy=mean(accuracyMeans),se=se(accuracyMeans))%>%ungroup()

`summarise()` regrouping output by 'nSolutions', 'pID' (override with `.groups` argument)
`summarise()` regrouping output by 'nSolutions' (override with `.groups` argument)

dataInput3$sol = as.factor(dataInput3$nSolutions>=1)
dataInput3$phaseT = recode(dataInput3$phaseT, '0' = "Low IC", '1' = "High IC")
dataInput3$nSolutions=dataInput3$nSolutions+1
  
plo=  ggplot(data=dataInput3, aes(y=accuracy, x=nSolutions)) +
  geom_errorbar(aes(ymin=accuracy-se, ymax=accuracy+se), width=.1)+
  geom_point(shape=23, size=3, fill="red")+
  labs(title="Accuracy and the Number of Solutions",
       x="log(Number of Solutions)",y="Accuracy")+
  coord_cartesian(ylim = c(0.3,1))+
  theme_light()+
  scale_x_log10()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5),strip.text.x = element_blank())+
  facet_grid(phaseT~sol, scales = "free_x", space = "free_x")+
  geom_smooth(data=dataInput3, aes(y=accuracy, x=nSolutions),method=glm,se=FALSE,fullrange=TRUE,linetype = "dashed")


outputName = 'tsp_nsol_acc_g'
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_nsol_acc_g"
quartz_off_screen 
                2 

dataInput=dataTrial_TSP

dataInput = dataInput %>% filter(sol==1)

dataInput$nSolutions_log =log(dataInput$nSolutions)
logitRandomIntercept = brm2(correct ~ nSolutions_log + phaseT + phaseT:nSolutions_log + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0)

tableName='tsp_acc_ns1'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

This calculates the significance of the p-value of the beta(slope) of number of witnesses for instances with high TCC.

fit = hypothesis(allModels[['tsp_acc_ns1']],"nSolutions_log + nSolutions_log:phaseT = 0", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |      95% HDI
------------------------
H1        | [0.31, 0.54]

Marginal Effects

plot(marginal_effects(allModels[['tsp_acc_ns1']]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

2.2.5.1 NSol (regressed alone)

#st02
dataInput=dataTrial_TSP
dataInput = dataInput %>% filter(sol==1)
dataInput$nSolutions_log =log(dataInput$nSolutions)

#Includes a dummy if nSolutions==0, that's variable: sol.
logitRandomIntercept = brm2(correct ~ nSolutions_log + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='tsp_acc_ns2'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

mean_acc = dataInput %>% group_by(instanceNumber,nSolutions_log, sol,phaseT) %>%
  summarise(accuracy=mean(correct)) 

`summarise()` regrouping output by 'instanceNumber', 'nSolutions_log', 'sol' (override with `.groups` argument)

logitRandomIntercept = allModels[['tsp_acc_ns2']]
pp  = plot(conditional_effects(logitRandomIntercept), plot = FALSE, ask = FALSE)

pp$nSolutions_log + 
  geom_point(data=mean_acc,aes(x = nSolutions_log, y = accuracy),inherit.aes = FALSE)

# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$nSolutions_log
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=mean_acc,
                               aes(x = nSolutions_log, y = accuracy, col=as.factor(sol),
                                   shape=as.factor(phaseT), size=2.5), inherit.aes = FALSE),
                         pp_plot$layers)


plo = pp_plot +
  xlab("Number of solution Witnesses (ln)")+#expression(IC[expost]))+
  ylab("Human Performance")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  scale_x_continuous(breaks = c(log(1),log(10),log(100),log(1000),log(10000)),
                     limits =c(0,NA),labels = c("1","10","100","1000","10000") )+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(0.2,1)


plo
ggsave(paste0(folderOut_figures,"/Nwit_TSP_acc.pdf"),plo,width = 6,height =6,units="in")

2.2.5.2 NSol (with TCC and no interaction)

#st03
dataInput=dataTrial_TSP
dataInput = dataInput %>% filter(sol==1)
dataInput$nSolutions_log =log(dataInput$nSolutions)


logitRandomIntercept = brm2(correct ~ nSolutions_log + phaseT + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName='tsp_acc_ns3'
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

2.2.6 IC_expost

dataInput = dataTrial_TSP

# Summarise the data to plot
mean_accuracy = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(accuracy = mean(correct))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

ggplot(mean_accuracy,aes(x = ICexpost, y = accuracy))+
  geom_point() +
  theme_light() +
  stat_smooth(formula =  y ~ I(x^0.01), method="lm", se= FALSE)+
  xlab("IC_expost")+
  ylab("Human Accuracy")

2.2.6.1 All instances

model_ICexpost = brm2(correct ~ ICexpost + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0)

tableName='tsp_acc_icexpost'
save_summarise_model(model_ICexpost, tableName)

Model Name already exists!

model_ICexpost = allModels[['tsp_acc_icexpost']]
pp  = plot(conditional_effects(model_ICexpost), plot = TRUE, ask = FALSE)

mean_accuracy$correct= mean_accuracy$accuracy
pp$ICexpost + 
  geom_point(data=mean_accuracy,aes(x = ICexpost, y = correct, col=as.factor(phaseT)),inherit.aes = FALSE)+
  theme_minimal()

# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$ICexpost
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=mean_accuracy,
                               aes(x = ICexpost, y = correct, col=as.factor(sol),
                                   shape=as.factor(phaseT), size=2.5), inherit.aes = FALSE),
                         pp_plot$layers)


plo = pp_plot +
  xlab("Instance Complexity (IC)")+#expression(IC[expost]))+
  ylab("Human Performance")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(0.2,1)

plo
ggsave(paste0(folderOut_figures,"/IC_TSP_acc.pdf"),plo,width = 6,height =6,units="in")

2.2.6.1.1 Model Fit

#st05
dataInput = dataTrial_TSP

# Make predictions excluding random effects (pID)
predictions = predict(model_ICexpost,dataInput, re_formula = NA)


# Estimating the significance of the fit. This is done considering the probability estimation rather than the binary classification.

#Performs the Hosmer-Lemeshow goodness of fit test
logistic_significance = generalhoslem::logitgof(dataInput$correct, predictions[,1], g = 10, ord = FALSE)

At least one cell in the expected frequencies table is < 1. Chi-square approximation may be incorrect.

logistic_significance

    Hosmer and Lemeshow test (binary model)

data:  dataInput$correct, predictions[, 1]
X-squared = 23.813, df = 8, p-value = 0.002463

# Finds R2 using binary outcomes
# https://stackoverflow.com/questions/40901445/function-to-calculate-r2-r-squared-in-r
#predictions = predict(model_ICexpost,dataInput, re_formula = NA)
r_2_binary = cor(dataInput$correct, predictions[,1])^2
r_2_binary

[1] 0.1659575

# Finds R2 using mean accuracies per instance
predictions2 = predict(model_ICexpost,mean_accuracy, re_formula = NA)
r_2_probabilities = cor(mean_accuracy$accuracy, predictions2[,1])^2
r_2_probabilities

[1] 0.7416523

2.2.6.2 Unsatisfiable instances

dataInput = dataTrial_TSP
dataInput = dataInput %>% filter(sol==0)
model_ICexpost = brm2(correct ~ ICexpost + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0)

tableName='tsp_acc_icexpost2'
save_summarise_model(model_ICexpost, tableName)

Model Name already exists!

#st05
dataInput = dataTrial_TSP %>% filter(sol==0)

model_ICexpost = allModels[['tsp_acc_icexpost2']]
# Make predictions excluding random effects (pID)
predictions = predict(model_ICexpost,dataInput, re_formula = NA)

# Estimating the significance of the fit. This is done considering the probability estimation rather than the binary classification.

#Performs the Hosmer-Lemeshow goodness of fit test
logistic_significance = generalhoslem::logitgof(dataInput$correct, predictions[,1], g = 10, ord = FALSE)

At least one cell in the expected frequencies table is < 1. Chi-square approximation may be incorrect.

logistic_significance

    Hosmer and Lemeshow test (binary model)

data:  dataInput$correct, predictions[, 1]
X-squared = 37.059, df = 8, p-value = 1.123e-05

# Finds R2 using binary outcomes
# https://stackoverflow.com/questions/40901445/function-to-calculate-r2-r-squared-in-r
#predictions = predict(model_ICexpost,dataInput, re_formula = NA)
r_2_binary = cor(dataInput$correct, predictions[,1])^2
r_2_binary

[1] 0.1382572

# Finds R2 using mean accuracies per instance

mean_accuracy_unsat = mean_accuracy %>% filter(sol==0)
predictions2 = predict(model_ICexpost,mean_accuracy_unsat, re_formula = NA)
r_2_probabilities = cor(mean_accuracy_unsat$accuracy, predictions2[,1])^2
r_2_probabilities

[1] 0.7899152

2.3 Time Spent

The default data used in this section is dataTrial_TSP_Time. This dataset excludes those participants that never skipped to answer submission.

We first calculate some summary stats for the whole data set (dataTrial_TSP) and see how many participants were excluded in dataTrial_TSP_Time.

2.3.0.1 SummStats TSP

#Summary Stats for Decision Problem
dataInput=dataTrial_TSP

timeSummary = dataInput %>% group_by(pID) %>% summarise(acc=mean(timeSpent)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() ungrouping output (override with .groups argument)

kable(timeSummary, digits=1, caption = "Time Summary")

Time Summary

mean	min	max	SD
32.2	19.9	39.2	5.2

yesNoProportions = dataInput %>% group_by(sol,pID) %>% summarise(acc=mean(timeSpent)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() regrouping output by ‘sol’ (override with .groups argument) summarise() ungrouping output (override with .groups argument)

kable(yesNoProportions, digits=1, caption = "Time Spent By Solution")

Time Spent By Solution

sol	mean	min	max	SD
0	34.1	20.0	40.0	5.9
1	30.2	19.9	38.5	5.2

hist(dataInput$timeSpent,breaks=40)

all_pIDs = unique(dataTrial_TSP$pID)
pIDs_not_in_Time = all_pIDs[!(all_pIDs %in% unique(dataTrial_TSP_Time$pID))]
print(paste0("The participants removed in the TIME df are ", length(pIDs_not_in_Time)))

[1] "The participants removed in the TIME df are 0"

obs_removed = nrow(dataTrial_TSP)-nrow(dataTrial_TSP_Time)
print(paste0("The number of observations removed in the TIME df is ", obs_removed))

[1] "The number of observations removed in the TIME df is 0"

2.3.0.2 SummStats TSP_Time

#Summary Stats for Decision Problem
dataInput=dataTrial_TSP_Time

timeSummary = dataInput %>% group_by(pID) %>% summarise(acc=mean(timeSpent)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() ungrouping output (override with .groups argument)

kable(timeSummary, digits=1, caption = "Time Summary")

Time Summary

mean	min	max	SD
32.2	19.9	39.2	5.2

yesNoProportions = dataInput %>% group_by(sol,pID) %>% summarise(acc=mean(timeSpent)) %>% summarise(mean=mean(acc),min=min(acc),max=max(acc),SD=sd(acc))

summarise() regrouping output by ‘sol’ (override with .groups argument) summarise() ungrouping output (override with .groups argument)

kable(yesNoProportions, digits=1, caption = "Time Spent By Solution")

Time Spent By Solution

sol	mean	min	max	SD
0	34.1	20.0	40.0	5.9
1	30.2	19.9	38.5	5.2

hist(dataInput$timeSpent,breaks=40)

2.3.0.3 Trial Number and Time Spent

#Trial (experience effect) effect on timeSpent
dataInput=dataTrial_TSP_Time
dataInput$totalTrial = dataInput$block * dataInput$trial

summaryByBlock = dataInput %>% group_by(block,pID) %>% summarise(time=mean(timeSpent)) %>% summarise(mean=mean(time),min=min(time),max=max(time),SD=sd(time))

`summarise()` regrouping output by 'block' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

kable(summaryByBlock,digits=2 , caption="Time Spent per Trial segregated By Block")

Time Spent per Trial segregated By Block

block	mean	min	max	SD
1	33.77	20.29	40.00	4.89
2	32.31	20.38	40.00	5.71
3	30.42	19.10	37.74	5.72

#Regression
linearRandomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ totalTrial+ (1|pID),
                            data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0, iter = iters_high)


tableName='tsp_time_01_r'
save_summarise_model(linearRandomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

2.3.1 TCC

dataInput=dataTrial_TSP_Time %>% group_by(phaseT,pID)%>%summarise(timeSpent1=mean(timeSpent))%>%
  ungroup()%>%group_by(phaseT)%>%summarise(timeSpent=mean(timeSpent1),se=se(timeSpent1))%>%ungroup()

`summarise()` regrouping output by 'phaseT' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

dataInput$phaseT = recode(dataInput$phaseT, '0' = "Low IC", '1' = "High IC")

plo= ggplot(data=dataInput, aes(y=timeSpent, x=as.factor(phaseT),label = round(timeSpent,digits = 1))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=timeSpent-se, ymax=timeSpent+se), width=.1)+
  #geom_signif(comparisons = list(c("Low IC","High IC")), annotations="***", y_position = 0.95, tip_length = 0.03)+
  labs(title="Time Spent In and Out of phase Transition",x="Instance complexity",y="Time Spent")+
  #coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white")+
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))

outputName = "tsp_time_02_g"
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_time_02_g"
quartz_off_screen 
                2 

#Anova contrast for In/Out Phase transition
dataInput = dataTrial_TSP_Time
anovaModel=aov(timeSpent~phaseT+Error(pID/phaseT),dataInput)
anovaPValue=summary(anovaModel)[[2]][[1]][['Pr(>F)']][[1]]
print(paste("P-value for one way ANOVA:",signif(anovaPValue,digits=3)))

[1] “P-value for one way ANOVA: 1.09e-05”

rm(dataInput)

timeSpent ~ phaseT + (1|pID)

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP_Time

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(timeSpent | cens(censored_time) ~ phaseT + (1|pID),
                       data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_time_02_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

2.3.1.0.1 Plot TCC

library(forcats)
library(ggpol)
#dataInput = dataTrial_TSP

# dataInput$timeSpent_pct = dataInput$timeSpent/SATTaskMaxTime
# dataInput2 = dataInput %>% 
#   group_by(instanceNumber,phaseT) %>% 
#   summarise(timeSpent_pct=median(timeSpent_pct))

dataInput = dataTrial_TSP_Time %>%
  mutate(region = fct_relevel(region,
                              "Underconstrained", "Phase Transition", "Overconstrained")) %>% 
  mutate(sol= as.factor(sol)) %>% 
  mutate(sol = fct_relevel(sol,"1", "0"))

dataInput$timeSpent_pct = dataInput$timeSpent/TSPTaskMaxTime
dataInput2 = dataInput %>% 
  group_by(instanceNumber,region,sol) %>% 
  summarise(timeSpent_pct=median(timeSpent_pct))

`summarise()` regrouping output by 'instanceNumber', 'region' (override with `.groups` argument)

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7


plo = ggplot(dataInput2,aes(x=region,y=timeSpent_pct,fill=sol))+
  # geom_boxjitter(aes(fill=region),jitter.color = NA,jitter.shape = 21)+
  geom_boxjitter(jitter.color = NA,jitter.shape = 21,
                 jitter.params = list(height=0,seed=10),
                 outlier.shape= NA)+
                 #,outlier.shape = 4, outlier.size=0.9)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  scale_fill_manual(name  ="Region",
                    breaks=c("1","0"),
                    labels=c("Satisfiable","Unsatisfiable"),
                    values=c("#90BE6D","#F94144"))+
  xlab("Typical Case Complexity (TCC)")+
  ylab("Time-on-task")+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(0.15,1)

plo
ggsave(paste0(folderOut_figures,"/TCC_TSP_time.pdf"),plo,width = 7,height =6,units="in")

2.3.2 Satisfiability and TCC

dataInput=dataTrial_TSP_Time %>% group_by(phaseT,sol,pID)%>%summarise(timeSpent1=mean(timeSpent))%>%ungroup()%>%group_by(sol,phaseT)%>%summarise(timeSpent=mean(timeSpent1),se=se(timeSpent1))%>%ungroup()

`summarise()` regrouping output by 'phaseT', 'sol' (override with `.groups` argument)
`summarise()` regrouping output by 'sol' (override with `.groups` argument)

sol_names <- c(
                    `0` = "Correct answer: NO",
                    `1` = "Correct answer: YES",
                    `2` = "If this is here it means data not filtered"
                    )

plo=ggplot(data=dataInput, aes(y=timeSpent, x=as.factor(as.logical(phaseT)), 
                               label = round(timeSpent,digits = 1),group=1)) +
  geom_line()+
  geom_errorbar(aes(ymin=timeSpent-se, ymax=timeSpent+se), width=.1) +
  geom_point(shape=21, size=3, fill="white")+
  labs(title="Time Spent segregated by solvanbility",x="In Phase Transition?",y="Time Spent")+
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle=0,hjust=0.5))+
  geom_text(hjust=-0.5, colour="black")+
  facet_grid(.~sol, labeller = as_labeller(sol_names))

outputName = "tsp_time_03_g"
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_time_03_g"
quartz_off_screen 
                2 

timeSpent ~ phaseT + sol + phaseT:sol + (1|pID)

dataInput= dataTrial_TSP_Time
#logistic with random effects (intercept)

dataInput$phaseT = as.factor(dataInput$phaseT)
dataInput$sol = as.factor(dataInput$sol)

# TTTODO
randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ phaseT + sol + phaseT:sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms,
                      refresh = 0, iter = iters_high)

tableName="tsp_time_03_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

Margianl Effects

#marginal_effects(randomIntercept)
plot(marginal_effects(allModels[["tsp_time_03_r"]]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

2.3.3 Satisfiability

#st06
dataInput= dataTrial_TSP_Time
#logistic with random effects (intercept)
#dataInput$sol =as.factor(dataInput$sol)

randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_time_03_r_B"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

2.3.4 Region effect on time-on-task (Constrainedness)

dataInput=dataTrial_TSP_Time 

dataInput2=dataInput %>% group_by(region,pID)%>%summarise(timeSpent1=mean(timeSpent))%>%ungroup()%>%group_by(region)%>%summarise(timeSpent=mean(timeSpent1),se=se(timeSpent1))%>%ungroup()

`summarise()` regrouping output by 'region' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

dataInput2$region <- factor(dataInput2$region, levels = c('Underconstrained',
                                                          'Phase Transition',
                                                          'Overconstrained'))

plo= ggplot(data=dataInput2, aes(y=timeSpent, x=region, label = round(timeSpent,digits = 1))) +
  geom_bar(stat = "identity")+
  geom_errorbar(aes(ymin=timeSpent-se, ymax=timeSpent+se), width=.1)+
  labs(title="time Spent by region",x="Region",y="time Spent")+
  geom_text(hjust=0.5,vjust = 5, colour="white") +
  theme_light()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))

outputName = "tsp_time_04_g"
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_time_04_g"
quartz_off_screen 
                2 

#Stats test for time in and out of phase transition: ANOVA
dataInput$overConstrained= (dataInput$region=='Overconstrained')
dataInput$underConstrained= (dataInput$region=='Underconstrained')

randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ overConstrained + underConstrained + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_time_04_g_A"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

Is there a difference between time-on-task between overconstrained and the underconstrained region?

fit = hypothesis(allModels[["tsp_time_04_g_A"]],"overConstrainedTRUE=underConstrainedTRUE", seed=seed_brms)
print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |      95% HDI
------------------------
H1        | [0.14, 0.21]

2.3.5 Number of solution witnesses that satisfy the xonstraints

dataInput=dataTrial_TSP_Time

#Plots nSolutions (x) vs. Accuracy (y)
dataInput3 = dataInput %>% group_by(nSolutions,pID,phaseT)%>%summarise(timeSpent1=mean(timeSpent))%>%ungroup()%>%group_by(nSolutions,phaseT)%>%summarise(timeSpent=mean(timeSpent1),se=se(timeSpent1))%>%ungroup()

`summarise()` regrouping output by 'nSolutions', 'pID' (override with `.groups` argument)
`summarise()` regrouping output by 'nSolutions' (override with `.groups` argument)

dataInput3$sol = dataInput3$nSolutions>=1
dataInput3$phaseT = recode(dataInput3$phaseT, '0' = "Low IC", '1' = "High IC")

dataInput3$nSolutions = dataInput3$nSolutions +1 

plo= ggplot(data=dataInput3, aes(y=timeSpent, x=nSolutions)) +
  geom_errorbar(aes(ymin=timeSpent-se, ymax=timeSpent+se), width=.1)+
  geom_point(shape=23, size=3, fill="red")+
  labs(title="Time Spent and the Number of Solutions",
       x="Number of   Solutions",y="Time Spent")+
  #coord_cartesian(ylim = c(0.5,1))+
  theme_light()+
  scale_x_log10()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5),strip.text.x = element_blank())+
  facet_grid(phaseT~sol, scales = "free_x", space = "free_x")+
  geom_smooth(data=dataInput3, aes(y=timeSpent, x=nSolutions),method=glm,
              se=FALSE,fullrange=TRUE,linetype = "dashed")

outputName = "tsp_nsols_time_g"
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_nsols_time_g"
quartz_off_screen 
                2 

dataInput=dataTrial_TSP_Time

dataInput = dataInput %>% filter(sol==1)
dataInput$nSolutions_log = log(dataInput$nSolutions)

linRandomIntercept = brm2(timeSpent_pct | cens(censored_time) ~ phaseT + nSolutions_log + phaseT:nSolutions_log + (1|pID),
                         data=dataInput,
                         chains=chains_brms,
                         cores = cores_brms,
                         seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_nsols_time_r"
save_summarise_model(linRandomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

This calculates the significance of the number of witnesses beta(slope) for instance with high TCC.

#This version recodes phaseT to OutphaseT, to 

fit = hypothesis(allModels[["tsp_nsols_time_r"]],"nSolutions_log + phaseT:nSolutions_log = 0", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

hdi(fit$samples,ci=0.95)

# Highest Density Interval

Parameter |        95% HDI
--------------------------
H1        | [-0.03, -0.01]

Marginal Effects

#marginal_effects(randomIntercept)
plot(marginal_effects(allModels[["tsp_nsols_time_r"]]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

2.3.5.1 No. of witnesses (Regressed alone)

#st07
dataInput= dataTrial_TSP_Time
dataInput = dataInput %>% filter(sol==1)
#logistic with random effects (intercept)
#dataInput$sol =as.factor(dataInput$sol)
dataInput$nSolutions_log = log(dataInput$nSolutions)

randomIntercept = brm2(timeSpent_pct | cens(censored_time) ~  nSolutions_log + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_nsols_time_r_B"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

#dataInput3
mean_acc = dataInput %>% group_by(instanceNumber,nSolutions,nSolutions_log, sol,phaseT) %>%
  summarise(accuracy=mean(correct),
            timeSpent_med =median(timeSpent)) 

`summarise()` regrouping output by 'instanceNumber', 'nSolutions', 'nSolutions_log', 'sol' (override with `.groups` argument)

logitRandomIntercept = allModels[['tsp_nsols_time_r_B']]
pp  = plot(conditional_effects(logitRandomIntercept), plot = FALSE, ask = FALSE)

pp$nSolutions + 
  geom_point(data=mean_acc,aes(x = nSolutions_log, y = timeSpent_med/TSPTaskMaxTime),inherit.aes = FALSE)

NA
NA

# Improving the plot

dataInput$nSolutions_log_bin = cut_width(dataInput$nSolutions_log, width=1, boundary =0,labels=FALSE)
dataInput$nSolutions_log_bin = (dataInput$nSolutions_log_bin*1)-0.5

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$nSolutions_log
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_boxplot(data=dataInput,
                                      aes(x = nSolutions_log_bin, 
                                          group = nSolutions_log_bin, 
                                          y = timeSpent_pct), 
                                      inherit.aes = FALSE,width=0.7),
                        geom_point(data=mean_acc, 
                                   aes(x = nSolutions_log, y = timeSpent_med/TSPTaskMaxTime,
                                       col=as.factor(sol), shape=as.factor(phaseT), size=2.5),
                                   inherit.aes = FALSE),
                         pp_plot$layers)

plo = pp_plot +
  xlab("Number of solution Witnesses")+#expression(IC[expost]))+
  ylab("Time-on-task")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+#c( "#FC4E07","#E7B800")
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  ylim(0, 1)+
  scale_x_continuous(breaks = c(log(1),log(10),log(100),log(1000),log(10000)),
                     limits =c(0,NA),labels = c("1","10","100","1000","10000") )+
  guides(shape = guide_legend(override.aes = list(size = 3)))

plo
ggsave(paste0(folderOut_figures,"/Nwit_TSP_time.pdf"),plo,width = 6,height =6,units="in")

# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$nSolutions_log
pp_plot$layers[[1]]$geom_params$se = TRUE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=dataInput,
                               aes(x = nSolutions_log, y = timeSpent_pct,
                                   shape=as.factor(phaseT), size=0.4,
                                   alpha = 0.5),
                                   inherit.aes = FALSE),
                    geom_jitter(data=mean_acc, 
                                   aes(x = nSolutions_log, 
                                       y = timeSpent_med/TSPTaskMaxTime,
                                       col=as.factor(sol),
                                       shape=as.factor(phaseT), size=2.5),
                                   inherit.aes = FALSE),
                         pp_plot$layers)

plo = pp_plot +
  xlab("Number of solution Witnesses")+#expression(IC[expost]))+
  ylab("Time-on-task")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+#c( "#FC4E07","#E7B800")
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  # ylim(0, 1.1)+
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),limits =c(0,1.1))+
  scale_x_continuous(breaks = c(log(1),log(10),log(100),log(1000),log(10000)),
                     limits =c(log(2),NA),labels = c("1","10","100","1000","10000") )+
  guides(shape = guide_legend(override.aes = list(size = 3)))

plo
ggsave(paste0(folderOut_figures,"/Nwit_TSP_time2.pdf"),plo,width = 6,height =6,units="in")

2.3.6 IC_expost

dataInput = dataTrial_TSP_Time

# Summarise the data to plot
mean_accuracy_time = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(accuracy = mean(correct), 
            timeSpent_med = median(timeSpent),
            timeSpent = mean(timeSpent))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

ggplot(mean_accuracy_time,aes(x = ICexpost, y = timeSpent))+
  geom_point() +
  theme_light() +
  stat_smooth(formula =  y ~ I(x^0.01), method="lm", se= FALSE)+
  xlab("IC_expost")+
  ylab("Time Spent")

model_ICexpost = brm2(timeSpent_pct | cens(censored_time) ~ ICexpost + (1|pID),
                           data=dataInput,
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms,
                           refresh = 0, iter = iters_high)

tableName='tsp_time_icexpost'
save_summarise_model(model_ICexpost, tableName)

Model Name already exists!Results may not be meaningful for censored models.Results may not be meaningful for censored models.Results may not be meaningful for censored models.

mean_accuracy_time$correct= mean_accuracy_time$accuracy

model_ICexpost = allModels[['tsp_time_icexpost']]
pp  = plot(conditional_effects(model_ICexpost), plot = FALSE, ask = FALSE)

dataInput$ICexpost_bin = cut_width(dataInput$ICexpost, width=0.02, boundary =0,labels=FALSE)
dataInput$ICexpost_bin = (dataInput$ICexpost_bin*0.02)-0.01

pp_plot = pp$ICexpost

pp_plot$layers <- c(geom_point(data=mean_accuracy_time,
                                   aes(x = ICexpost, y = timeSpent_med/TSPTaskMaxTime),
                                   inherit.aes = FALSE),
                         pp_plot$layers)

pp_plot

pp_plot = pp$ICexpost

pp_plot$layers[[1]]$geom_params$se = TRUE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_point(data=dataInput,
                               aes(x = ICexpost, y = timeSpent_pct,
                                   shape=as.factor(phaseT), size=0.4,
                                   #col=as.factor(sol),
                                   alpha = 0.3),
                                   inherit.aes = FALSE),
                    geom_jitter(data=mean_accuracy_time, 
                                   aes(x = ICexpost, 
                                       y = timeSpent_med/TSPTaskMaxTime,
                                       col=as.factor(sol),
                                       shape=as.factor(phaseT), size=2.5),
                                   inherit.aes = FALSE),
                         pp_plot$layers)

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

plo = pp_plot +
  xlab("Instance Complexity (IC)")+#expression(IC[expost]))+
  ylab("Time-on-task")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+#c( "#FC4E07","#E7B800")
  # scale_shape_manual(name="IC",values = c(2, 8))+
  # scale_color_manual(name="Solution",values = c("red", "blue"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1),limits =c(0,1.1))+
  guides(shape = guide_legend(override.aes = list(size = 3)))

plo
ggsave(paste0(folderOut_figures,"/IC_TSP_time1.pdf"),plo,width = 6,height =6,units="in")

2.4 Time-on-task and accuracy

dataInput= dataTrial_TSP_Time %>% group_by(correct,pID) %>% summarise(timeSpentMeans=mean(timeSpent)) %>% ungroup() %>% group_by(correct) %>%
  summarise(time=mean(timeSpentMeans),se=se(timeSpentMeans))%>%ungroup()

`summarise()` regrouping output by 'correct' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

plo = ggplot(data=dataInput, aes(y=time, x=as.factor(as.logical(correct)), label = round(time,digits = 2),group=1)) +
  geom_line()+
  geom_errorbar(aes(ymin=time-se, ymax=time+se), width=.1) +
  geom_point(shape=21, size=3, fill="white")+
  labs(title="Time Spent on Correct/Incorrect instances",x="Answer Correct?",y="Time Spent") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle=0,hjust=0.5))+
  geom_text(hjust=-0.5, colour="black")

outputName = "tsp_time_08_g_B"
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_time_08_g_B"
quartz_off_screen 
                2 

correct ~ timeSpent + (1|pID)

dataInput = dataTrial_TSP_Time
#Stats test for time in and out of phase transition: ANOVA
anovaModel=aov(timeSpent~correct+Error(pID/correct),dataInput)
anovaPValue=summary(anovaModel)[[2]][[1]][['Pr(>F)']][[1]]
print(paste("P-value for one way ANOVA:",signif(anovaPValue,digits=3)))

[1] "P-value for one way ANOVA: 4.43e-08"

#Stats test for time for correct/incorrect answers
logitRandomIntercept = brm2(correct ~ timeSpent_pct + (1|pID),
                           family=bernoulli(link="logit"),
                           data=dataInput, 
                           chains=chains_brms,
                           cores = cores_brms,
                           seed=seed_brms, refresh = 0)

tableName="tsp_time_08_r_B"
save_summarise_model(logitRandomIntercept, tableName)

Model Name already exists!

Marginal Effects {#tsp_time_acc}

plot(marginal_effects(allModels[["tsp_time_08_r_B"]]), plot = TRUE, ask = FALSE)

Method 'marginal_effects' is deprecated. Please use 'conditional_effects' instead.

2.5 Number of Clicks

2.5.0.1 Satisfiability

dataInput=dataTrial_TSP_clicks

# BOXPLOT
plo= ggplot(data=dataInput, aes(y=n_clicks, x=region, label = round(n_clicks,digits = 0))) +
  geom_boxplot()+
  labs(x="Region",y="n_clicks")+
  #coord_cartesian(ylim = c(0.5,1))+
  geom_text(hjust=0.5,vjust = 5, colour="white") +
  theme_light()+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),plot.title = element_text(hjust = 0.5))+
  facet_grid(.~sol)


outputName = 'tsp_clicks_00_g'
plotExport(plo,outputName,folderOut_figures)

[1] "tsp_clicks_00_g"
quartz_off_screen 
                2 

size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

dataInput2 = dataTrial_TSP_clicks %>% 
  mutate(region = fct_relevel(region,
                              "Underconstrained", "Phase Transition", "Overconstrained"))%>% 
  mutate(sol= as.factor(sol)) %>% 
  mutate(sol = fct_relevel(sol,"1", "0"))

dataInput3 = dataInput2 %>% 
  group_by(instanceNumber,region,sol) %>% 
  summarise(n_clicks=mean(n_clicks))

`summarise()` regrouping output by 'instanceNumber', 'region' (override with `.groups` argument)

plo = ggplot(dataInput3,aes(x=region,y=n_clicks,fill=sol))+
  # geom_boxjitter(aes(fill=region),jitter.color = NA,jitter.shape = 21)+
  geom_boxjitter(jitter.color = NA,jitter.shape = 21,
                 jitter.params = list(height=0,seed=10),
                 outlier.shape= NA)+
                 #,outlier.shape = 4, outlier.size=0.9)+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5))+
  scale_fill_manual(name  ="Region",
                    breaks=c("1","0"),
                    labels=c("Satisfiable","Unsatisfiable"),
                    values=c("#90BE6D","#F94144"))+
  xlab("Typical Case Complexity (TCC)")+
  ylab("Number of clicks")+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  ylim(20,30)


plo
ggsave(paste0(folderOut_figures,"/TCC_TSP_nclicks.pdf"),plo,width = 7,height =6,units="in")

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks ~ sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_clicks_01_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!

2.5.0.2 TCC

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks ~ phaseT + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_clicks_02_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!

2.5.0.2.1 And satisfiability

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks ~ region+ (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_clicks_02B_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP_clicks

#logistic with random effects (intercept)
# TTTODO
dataInput$phaseT = as.factor(dataInput$phaseT)
dataInput$sol = as.factor(dataInput$sol)
randomIntercept = brm2(n_clicks ~ phaseT + sol + phaseT:sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_clicks_02C_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!

pp  = plot(conditional_effects(allModels[["tsp_clicks_02C_r"]]), plot = TRUE, ask = FALSE)

# stche
model = allModels[['tsp_clicks_02C_r']]

fit = hypothesis(model,"sol1 + phaseT1:sol1=0", seed=seed_brms)

print(fit)

Hypothesis Tests for class b:

---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

2.5.0.3 ICexpost

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks ~ ICexpost + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_clicks_03_r"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!

pp  = plot(conditional_effects(allModels[["tsp_clicks_03_r"]]), plot = TRUE, ask = FALSE)

dataInput= dataTrial_TSP_clicks
# Summarise the data to plot
mean_nclicks = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions) %>% 
  summarise(n_clicks = median(n_clicks))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

model_ICexpost = allModels[["tsp_clicks_03_r"]]
pp  = plot(conditional_effects(model_ICexpost), plot = FALSE, ask = FALSE)

pp$ICexpost +
  geom_point(data=mean_nclicks,aes(x = ICexpost, y =n_clicks),inherit.aes = FALSE)

# Improving the plot
# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$ICexpost
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"



pp_plot$layers <- c(geom_jitter(data=dataInput,
                               aes(x = ICexpost, y = n_clicks,
                                   shape=as.factor(phaseT)), size=0.7,
                                   alpha = 0.3,
                                   inherit.aes = FALSE,height=0.2),
                    geom_point(data=mean_nclicks, 
                                   aes(x = ICexpost, y = n_clicks,
                                       col=as.factor(sol), shape=as.factor(phaseT)),size=2.5,
                                   inherit.aes = FALSE),
                         pp_plot$layers)

plo = pp_plot +
  xlab("Instance Complexity (IC)")+#expression(IC[expost]))+
  ylab("Number of Clicks")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  ylim(10,50)
#+  guides(shape = guide_legend(override.aes = list(size = 3)))

# Justification to reduce range to 10-50
sum(between(dataTrial_TSP_clicks$n_clicks,10,50))/length(dataTrial_TSP_clicks$n_clicks)

[1] 0.9873016

plo
ggsave(paste0(folderOut_figures,"/IC_TSP_nclicks.pdf"),plo,width = 6,height =6,units="in")

2.5.0.3.1 And satisfiability

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks ~ ICexpost + sol + ICexpost:sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_clicks_03_rB"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!

#Logistic Regression for In/Out Phase transition
dataInput= dataTrial_TSP_clicks

#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks ~ ICexpost + sol + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_clicks_03_rC"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!

pp  = plot(conditional_effects(allModels[["tsp_clicks_03_rC"]]), plot = TRUE, ask = FALSE)

2.5.0.4 Number of witnesses

Only Satisfiable instances:

dataInput= dataTrial_TSP_clicks %>% filter(sol==1)

dataInput$nSolutions_log = log(dataInput$nSolutions)
#logistic with random effects (intercept)
# TTTODO
randomIntercept = brm2(n_clicks ~ nSolutions_log + (1|pID),
                      data=dataInput, 
                      chains=chains_brms,
                      cores = cores_brms,
                      seed=seed_brms, refresh = 0, iter = iters_high)

tableName="tsp_clicks_04_rB"
save_summarise_model(randomIntercept, tableName)

Model Name already exists!

#dataInput = dataTrial_SAT_clicks
# Summarise the data to plot
mean_nclicks = dataInput %>% 
  group_by(instanceNumber, sol, ICexpost, phaseT, nSolutions_log) %>% 
  summarise(n_clicks = median(n_clicks))

`summarise()` regrouping output by 'instanceNumber', 'sol', 'ICexpost', 'phaseT' (override with `.groups` argument)

model_nsol = allModels[["tsp_clicks_04_rB"]]
pp  = plot(conditional_effects(model_nsol), plot = FALSE, ask = FALSE)

pp$nSolutions_log +
   geom_point(data=mean_nclicks, aes(x = nSolutions_log, y =n_clicks),inherit.aes = FALSE)

# Improving the plot
# Improving the plot
size_big = 20
size_small = 16
size_ss = 10
size_xs = 7

pp_plot = pp$nSolutions
pp_plot$layers[[1]]$geom_params$se = FALSE
pp_plot$layers[[1]]$aes_params$colour="#577590"

pp_plot$layers <- c(geom_jitter(data=dataInput,
                               aes(x = nSolutions_log, y = n_clicks,
                                   shape=as.factor(phaseT)), size=1,
                                   alpha = 0.3,
                                   inherit.aes = FALSE,width=0.2,height =0.2),
                        geom_point(data=mean_nclicks, 
                                   aes(x = nSolutions_log, y = n_clicks,
                                       col=as.factor(sol), shape=as.factor(phaseT)), size=4,
                                   inherit.aes = FALSE),
                         pp_plot$layers)

plo = pp_plot +
  xlab("Number of Solution Witnesses")+
  ylab("Number of Clicks")+
  scale_shape_manual(name="",values = c(17,16)) +
  scale_color_manual(name="",values = c("1"="#90BE6D","0"="#F94144"))+
  theme_classic()+
  theme(axis.title = element_text(size= size_big),
       axis.text=element_text(size=size_small),
       legend.position = "none")+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  scale_x_continuous(breaks = c(log(1),log(10),log(100),log(1000),log(10000)),
                   limits =c(0,NA),labels = c("1","10","100","1000","10000") )+
  ylim(10,50)

plo

ggsave(paste0(folderOut_figures,"/Nwit_TSP_nclicks.pdf"),
       plo,width = 6,height =6,units="in")

3 Save Models

#saveRDS(allModels,file.path(folderOut,"ICx3_models.RData"))

4 Session Info

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggpol_0.0.7        forcats_0.5.0      performance_0.4.8  bayestestR_0.8.2  
 [5] see_0.6.2          parameters_0.12.0  tidyr_1.1.1        DiagrammeR_1.0.6.1
 [9] bmlm_1.3.11        bayesplot_1.7.2    sjstats_0.18.0     brms_2.13.5       
[13] Rcpp_1.0.5.2       readr_1.3.1        Hmisc_4.4-1        Formula_1.2-3     
[17] survival_3.1-12    lattice_0.20-41    plotly_4.9.2.1     pander_0.6.3      
[21] ggsignif_0.6.0     knitr_1.29         ggplot2_3.3.2      reshape2_1.4.4    
[25] dplyr_1.0.2       

loaded via a namespace (and not attached):
  [1] backports_1.1.9      plyr_1.8.6           igraph_1.2.5        
  [4] lazyeval_0.2.2       splines_4.0.2        crosstalk_1.1.0.1   
  [7] rstantools_2.1.1     inline_0.3.15        digest_0.6.25       
 [10] htmltools_0.5.0      rsconnect_0.8.16     fansi_0.4.1         
 [13] magrittr_1.5         checkmate_2.0.0      cluster_2.1.0       
 [16] modelr_0.1.8         RcppParallel_5.0.2   matrixStats_0.56.0  
 [19] xts_0.12-0           prettyunits_1.1.1    jpeg_0.1-8.1        
 [22] colorspace_1.4-1     xfun_0.16            callr_3.4.3         
 [25] crayon_1.3.4         jsonlite_1.7.0       lme4_1.1-23         
 [28] zoo_1.8-8            glue_1.4.1           gtable_0.3.0        
 [31] emmeans_1.5.0        sjmisc_2.8.5         V8_3.2.0            
 [34] pkgbuild_1.1.0       rstan_2.21.2         abind_1.4-5         
 [37] scales_1.1.1         mvtnorm_1.1-1        miniUI_0.1.1.1      
 [40] viridisLite_0.3.0    xtable_1.8-4         generalhoslem_1.3.4 
 [43] htmlTable_2.0.1      tmvnsim_1.0-2        foreign_0.8-80      
 [46] stats4_4.0.2         StanHeaders_2.21.0-6 DT_0.15             
 [49] htmlwidgets_1.5.1    httr_1.4.2           threejs_0.3.3       
 [52] lavaan_0.6-8         RColorBrewer_1.1-2   ellipsis_0.3.1      
 [55] farver_2.0.3         pkgconfig_2.0.3      reshape_0.8.8       
 [58] loo_2.3.1            nnet_7.3-14          labeling_0.3        
 [61] tidyselect_1.1.0     rlang_0.4.7          later_1.1.0.1       
 [64] effectsize_0.3.2     visNetwork_2.0.9     munsell_0.5.0       
 [67] tools_4.0.2          cli_2.0.2            generics_0.0.2      
 [70] sjlabelled_1.1.6     broom_0.7.0          ggridges_0.5.2      
 [73] fdrtool_1.2.16       evaluate_0.14        stringr_1.4.0       
 [76] fastmap_1.0.1        yaml_2.2.1           processx_3.4.3      
 [79] purrr_0.3.4          glasso_1.11          pbapply_1.4-3       
 [82] nlme_3.1-148         mime_0.9             compiler_4.0.2      
 [85] shinythemes_1.1.2    rstudioapi_0.11      curl_4.3            
 [88] png_0.1-7            tibble_3.0.3         statmod_1.4.34      
 [91] pbivnorm_0.6.0       stringi_1.4.6        highr_0.8           
 [94] ps_1.3.4             qgraph_1.6.9         Brobdingnag_1.2-6   
 [97] Matrix_1.2-18        psych_2.0.9          nloptr_1.2.2.2      
[100] markdown_1.1         shinyjs_1.1          vctrs_0.3.2         
[103] pillar_1.4.6         lifecycle_0.2.0      bridgesampling_1.0-0
[106] estimability_1.3     corpcor_1.6.9        data.table_1.13.0   
[109] insight_0.13.1       httpuv_1.5.4         R6_2.4.1            
[112] latticeExtra_0.6-29  promises_1.1.1       gridExtra_2.3       
[115] codetools_0.2-16     boot_1.3-25          colourpicker_1.0    
[118] MASS_7.3-51.6        gtools_3.8.2         assertthat_0.2.1    
[121] withr_2.2.0          mnormt_2.0.2         shinystan_2.5.0     
[124] mgcv_1.8-31          hms_0.5.3            grid_4.0.2          
[127] rpart_4.1-15         coda_0.19-3          minqa_1.2.4         
[130] rmarkdown_2.7        lubridate_1.7.9      shiny_1.5.0         
[133] base64enc_0.1-3      dygraphs_1.1.1.6